diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2ee0be1c437..4ead8403e2f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/NonlinearSolverEbos.cpp opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp opm/simulators/flow/ValidationFunctions.cpp + opm/simulators/flow/aspinPartition.cpp opm/simulators/linalg/bda/WellContributions.cpp opm/simulators/linalg/bda/MultisegmentWellContribution.cpp opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp @@ -226,6 +227,7 @@ endif() # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES tests/test_ALQState.cpp + tests/test_aspinPartition.cpp tests/test_blackoil_amg.cpp tests/test_convergenceoutputconfiguration.cpp tests/test_convergencereport.cpp @@ -343,6 +345,7 @@ list (APPEND TEST_DATA_FILES tests/include/summary.inc tests/include/test1_20x30x10.grdecl tests/include/well_vfp.ecl + tests/partition.txt ) @@ -362,6 +365,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp opm/simulators/flow/KeywordValidation.hpp opm/simulators/flow/ValidationFunctions.hpp + opm/simulators/flow/aspinPartition.hpp + opm/simulators/flow/SubDomain.hpp opm/core/props/BlackoilPhases.hpp opm/core/props/phaseUsageFromDeck.hpp opm/core/props/satfunc/RelpermDiagnostics.hpp @@ -421,6 +426,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/SmallDenseMatrixUtils.hpp opm/simulators/linalg/WellOperators.hpp opm/simulators/linalg/WriteSystemMatrixHelper.hpp + opm/simulators/linalg/extractMatrix.hpp opm/simulators/linalg/findOverlapRowsAndColumns.hpp opm/simulators/linalg/getQuasiImpesWeights.hpp opm/simulators/linalg/setupPropertyTree.hpp diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 51d50cb6113..566cb851910 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -389,14 +389,16 @@ public: // of threshold pressure is a quite big hack that only makes sense for ECL // datasets. (and even there, its physical justification is quite // questionable IMO.) - if (std::abs(Toolbox::value(pressureDifference)) > thpres) { - if (pressureDifference < 0.0) - pressureDifference += thpres; - else - pressureDifference -= thpres; - } - else { - pressureDifference = 0.0; + if (thpres > 0.0) { + if (std::abs(Toolbox::value(pressureDifference)) > thpres) { + if (pressureDifference < 0.0) + pressureDifference += thpres; + else + pressureDifference -= thpres; + } + else { + pressureDifference = 0.0; + } } } diff --git a/ebos/eclgenericcpgridvanguard.cc b/ebos/eclgenericcpgridvanguard.cc index 2fc206bf5ac..3e23eb41724 100644 --- a/ebos/eclgenericcpgridvanguard.cc +++ b/ebos/eclgenericcpgridvanguard.cc @@ -49,7 +49,7 @@ #include #include #include -#endif //HAVE_DUNE_FEM +#endif // HAVE_DUNE_FEM #include #include @@ -84,6 +84,20 @@ void EclGenericCpGridVanguard::releaseEquilGrid() } #if HAVE_MPI + +template +std::vector +EclGenericCpGridVanguard:: +partitionCells(const int numDomains, + const double imbalanceTolerance, + const std::vector& wells) const +{ + const auto* faceTrans = static_cast(nullptr); + + return this->grid_-> + zoltanPartitionWithoutScatter(&wells, faceTrans, numDomains, imbalanceTolerance); +} + template void EclGenericCpGridVanguard:: doLoadBalance_(const Dune::EdgeWeightMethod edgeWeightsMethod, diff --git a/ebos/eclgenericcpgridvanguard.hh b/ebos/eclgenericcpgridvanguard.hh index c873b410b83..fa9f98ce598 100644 --- a/ebos/eclgenericcpgridvanguard.hh +++ b/ebos/eclgenericcpgridvanguard.hh @@ -28,9 +28,12 @@ #define EWOMS_ECL_CP_GRID_GENERIC_VANGUARD_HH #include + #include +#include #include +#include #include #include @@ -115,6 +118,13 @@ public: return this->cell_part_; } +#if HAVE_MPI + std::vector + partitionCells(const int numDomains, + const double imbalanceTolerance, + const std::vector& wells) const; +#endif + protected: /*! * \brief Distribute the simulation grid over multiple processes diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index e41959d0068..9aea91b8dac 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -1698,6 +1698,9 @@ public: */ bool recycleFirstIterationStorage() const { + if (this->simulator().timeStepIndex() == 0) + return true; + int episodeIdx = this->episodeIndex(); return !this->drsdtActive_(episodeIdx) && !this->drvdtActive_(episodeIdx) && diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index fcd6bd1edf2..c32f6ea393f 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -27,6 +27,24 @@ #include #include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + #include #include #include @@ -43,6 +61,9 @@ #include #include #include + + + #include #include #include @@ -52,17 +73,24 @@ #include #include + +#include +#include + #include #include #include #include +#include + #include #include #include #include #include #include +#include #include #include @@ -200,6 +228,8 @@ namespace Opm { typedef typename SparseMatrixAdapter::IstlMatrix Mat; typedef Dune::BlockVector BVector; + using Domain = SubDomain; + typedef ISTLSolverEbos ISTLSolverType; class ComponentName @@ -295,14 +325,104 @@ namespace Opm { // compute global sum of number of cells global_nc_ = detail::countGlobalCells(grid_); convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves. + // TODO: remember to fix! + if (param_.nonlinear_solver_ == "aspin" || param_.nonlinear_solver_ == "nldd" || param_.nonlinear_solver_ == "purelocal") { + setupAspinDomains(); + } } + bool isParallel() const { return grid_.comm().size() > 1; } + const EclipseState& eclState() const { return ebosSimulator_.vanguard().eclState(); } + + + void setupAspinDomains() + { + // Create partitions. + const auto& [partition_vector, num_domains] = + this->partitionCells(); + + // 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(cell == num_cells); + assert(count == sizes); + + // Create the domains. + //const int num_domains = partitions.size(); + 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.json"; + } 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) @@ -335,6 +455,10 @@ namespace Opm { std::cout << "equation scaling not supported yet" << std::endl; //updateEquationsScaling(); } + + // Setup domain->well mapping. + wellModel().setupDomains(domains_); + report.pre_post_time += perfTimer.stop(); return report; @@ -348,18 +472,11 @@ namespace Opm { /// \param[in] iteration should be 0 for the first call of a new timestep /// \param[in] timer simulation timer /// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control) - /// \param[in, out] reservoir_state reservoir state variables - /// \param[in, out] well_state well state variables template SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface& timer, NonlinearSolverType& nonlinear_solver) { - SimulatorReportSingle report; - failureReport_ = SimulatorReportSingle(); - Dune::Timer perfTimer; - - perfTimer.start(); if (iteration == 0) { // For each iteration we store in a vector the norms of the residual of // the mass balance for each active phase, the well flux and the well equations. @@ -370,8 +487,32 @@ namespace Opm { convergence_reports_.back().report.reserve(11); } + if (iteration == 0) { + return nonlinearIterationNewton(iteration, timer, nonlinear_solver); + } + if (param_.nonlinear_solver_ == "aspin" || param_.nonlinear_solver_ == "nldd" || param_.nonlinear_solver_ == "purelocal") { + return nonlinearIterationAspin(iteration, timer, nonlinear_solver); + } else { + return nonlinearIterationNewton(iteration, timer, nonlinear_solver); + } + } + + + template + SimulatorReportSingle nonlinearIterationNewton(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + + // ----------- Set up reports and timer ----------- + SimulatorReportSingle report; + failureReport_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + + perfTimer.start(); report.total_linearizations = 1; + // ----------- Assemble ----------- try { report += assembleReservoir(timer, iteration); report.assemble_time += perfTimer.stop(); @@ -383,12 +524,41 @@ namespace Opm { throw; // continue throwing the stick } +#if 0 + // DEBUG HACK + // Check if domain-diagonal parts are identical to stored subdomain matrices. + if (param_.nonlinear_solver_ == "nldd") { + const Mat& main_matrix = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + for (const auto& domain : domains_) { + if (domain_matrices_[domain.index]) { + auto domain_diag_part = Details::extractMatrix(main_matrix, domain.cells); + if (Details::matrixEqual(domain_diag_part, *domain_matrices_[domain.index])) { + OpmLog::debug(" *** matrices IDENTICAL for domain " + std::to_string(domain.index)); + } else { + OpmLog::debug(" *** matrices DIFFER for domain " + std::to_string(domain.index)); + } + } else { + OpmLog::debug(" *** no local matrix for domain " + std::to_string(domain.index)); + } + } + } +#endif + + // ----------- 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); + auto convrep = getConvergence(timer, iteration, residual_norms); + // if (!convrep.converged()) { + // std::ostringstream os; + // os << "Convergence data for Newton iteration " << iteration << "\n" + // << convrep; + // OpmLog::debug(os.str()); + // } + ////wellModel().logPrimaryVars(); + report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();; ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); convergence_reports_.back().report.push_back(std::move(convrep)); @@ -400,8 +570,11 @@ namespace Opm { OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); } } - report.update_time += perfTimer.stop(); + //report.convergence_check_time += perfTimer.stop(); residual_norms_history_.push_back(residual_norms); + // Dune::writeMatrixToMatlab(ebosSimulator_.model().linearizer().jacobian().istlMatrix(), "beforeWellModelLinearize.matlab"); + + // ----------- If not converged, solve linear system ----------- if (!report.converged) { perfTimer.reset(); perfTimer.start(); @@ -422,6 +595,7 @@ namespace Opm { // Note that linearize may throw for MSwells. wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), ebosSimulator().model().linearizer().residual()); + //wellModel().logPrimaryVars(); solveJacobianSystem(x); report.linear_solve_setup_time += linear_solve_setup_time_; @@ -444,12 +618,13 @@ namespace Opm { // on observation to avoid some big performance degeneration under some circumstances. // there is no theorectical explanation which way is better for sure. wellModel().postSolve(x); + //wellModel().logPrimaryVars(); if (param_.use_update_stabilization_) { // Stabilize the nonlinear update. bool isOscillate = false; bool isStagnate = false; - nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate); if (isOscillate) { current_relaxation_ -= nonlinear_solver.relaxIncrement(); current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); @@ -472,6 +647,731 @@ namespace Opm { return report; } + + + + template + SimulatorReportSingle nonlinearIterationAspin(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + // ----------- Set up reports and timer ----------- + SimulatorReportSingle report; + failureReport_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + + //wellModel().logPrimaryVars(); + + 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; + // todo (?): make the report an attribute of the class + throw; // continue throwing the stick + } + //wellModel().logPrimaryVars(); + + // ----------- 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); + if (!convrep.converged()) { + // std::ostringstream os; + // os << "Convergence data for first global iteration:\n" + // << convrep; + // OpmLog::debug(os.str()); + } + report.converged = convrep.converged() && iteration > nonlinear_solver.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.convergence_check_time += perfTimer.stop(); + residual_norms_history_.push_back(residual_norms); + + if (report.converged) { + return report; + } + + // ----------- If not converged, do an ASPIN or NLDD iteration ----------- + + // apply the Schur compliment of the well model to the reservoir linearized + // equations + // wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), + // ebosSimulator().model().linearizer().residual()); + + // Take a copy of the FI residual. + // auto fully_implicit_residual = ebosSimulator().model().linearizer().residual(); + // ... and the Jacobian. + // auto fully_implicit_jacobian = ebosSimulator().model().linearizer().jacobian().istlMatrix(); + // ... and the solution state that generated it. + auto& solution = ebosSimulator().model().solution(0); + auto initial_solution = solution; + // auto initial_wellstate = wellModel().wellState(); + auto locally_solved = initial_solution; + + // ----------- Decide on an ordering for the domains ----------- + std::vector domain_order(domains_.size()); + if (param_.local_solve_approach_ == "gauss-seidel") { + if (true) { + // 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; + } + } else { + // 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; + } + } + } else { + std::iota(domain_order.begin(), domain_order.end(), 0); + } + + // ----------- 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; + if (param_.local_solve_approach_ == "jacobi") { + auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; +#if 0 + 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 + if (/*!local_report.converged*/ false) { + // Try again with a less strict tolerance. + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + // auto param_orig = param_; + // param_.local_tolerance_scaling_cnv_ *= 10.0; + // param_.local_tolerance_scaling_mb_ *= 10.0; + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + // param_ = param_orig; + } + 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); + } +#endif + } else { + assert(param_.local_solve_approach_ == "gauss-seidel"); + auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + if (/*!local_report.converged*/ false) { + // Try again with a less strict tolerance. + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + // auto param_orig = param_; + // param_.local_tolerance_scaling_cnv_ = 50.0; + // param_.local_tolerance_scaling_mb_ = 50.0; + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + // param_ = param_orig; + } + 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& rf : convrep.reservoirFailures()) { + if (rf.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { + // mb_sum += rf.magnitude(); + } else if (rf.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { + // cnv_sum += rf.magnitude(); + } + } + // If not too high, we overrule the convergence failure. + if (mb_sum < 1e-3 && cnv_sum < 1.0) { + 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); + } + } + // 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; + } + +#if 0 + // DEBUG write solution vectors. + { + const int nc = UgGridHelpers::numCells(grid_); + BVector sol(nc); + for (int c = 0; c < nc; ++c) { + sol[c] = solution[c]; + } + BVector loc_sol(nc); + for (int c = 0; c < nc; ++c) { + loc_sol[c] = locally_solved[c]; + } + Dune::storeMatrixMarket(sol, param_.local_solve_approach_ + "-" + std::to_string(iteration) + ".mm"); + Dune::storeMatrixMarket(loc_sol, param_.local_solve_approach_ + "-local-" + std::to_string(iteration) + ".mm"); + } +#endif + + // Print summary of local solve convergence. + { + 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; + } + +#define ASPIN_EXTRA_WELL_OUTPUT 1 +#if ASPIN_EXTRA_WELL_OUTPUT + // Extra debug output. + //wellModel().logPrimaryVars(); + { + std::ostringstream os; + os << " ** Well rates:"; + for (int w = 0; w < wellModel().wellState().numWells(); ++w) { + os << " | "; + const auto& wr = wellModel().wellState().wellRates(w); + for (const double r : wr) { + os << ' ' << r; + } + } + OpmLog::debug(os.str()); + } +#endif + + if (param_.nonlinear_solver_ == "purelocal") { + if (param_.local_solve_approach_ == "jacobi") { + solution = locally_solved; + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + return report; + } + else if (param_.nonlinear_solver_ == "nldd") { + if (param_.local_solve_approach_ == "jacobi") { + solution = locally_solved; + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + // Finish with a Newton step. + // ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + auto rep = nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver); + report += rep; + if (rep.converged) { + report.converged = true; + } + return report; + } + report.total_newton_iterations = 1; + + /* + // HACK to check FI convergence + // Take a copy of the FI residual. + auto res1 = ebosSimulator().model().linearizer().residual(); + // ... and the Jacobian. + auto jac1 = ebosSimulator().model().linearizer().jacobian().istlMatrix(); + // ... and the solution state that generated it. + auto sol1 = solution; + auto ws1 = wellModel().wellState(); + + // ----------- Assemble ----------- + try { + report += assembleReservoir(timer, iteration); + report.assemble_time += perfTimer.stop(); + } + catch (...) { + report.assemble_time += perfTimer.stop(); + failureReport_ += report; + // todo (?): make the report an attribute of the class + throw; // continue throwing the stick + } + //wellModel().logPrimaryVars(); + + // ----------- Check if converged ----------- + 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 > nonlinear_solver.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(NumericalProblem, "Too large residual found!"); + } + } + report.update_time += perfTimer.stop(); + residual_norms_history_.push_back(residual_norms); + + if (report.converged) { + return report; + } + // HACK reinstate everything + ebosSimulator().model().linearizer().residual() = res1; + ebosSimulator().model().linearizer().jacobian().istlMatrix() = jac1; + solution = sol1; + wellModel().wellState() = ws1; + */ + + // ----------- Compute ASPIN residual, check convergence ----------- + const int nc = grid_.size(0); // Includes both owned cells and overlap. + BVector aspin_residual(nc); + const int num_vars = aspin_residual[0].size(); + std::vector max_sol(num_vars, 0.0); + std::vector max_diff(num_vars, 0.0); + for (int c = 0; c < nc; ++c) { + aspin_residual[c] = initial_solution[c] - locally_solved[c]; + for (int var = 0; var < num_vars; ++var) { + max_sol[var] = std::max(max_sol[var], std::fabs(locally_solved[c][var])); + max_diff[var] = std::max(max_diff[var], std::fabs(aspin_residual[c][var])); + } + } + double max_scaled_diff = 0.0; + for (int var = 0; var < num_vars; ++var) { + max_scaled_diff = std::max(max_scaled_diff, max_diff[var]/max_sol[var]); + } + OpmLog::debug(" Scaled ASPIN residual: " + std::to_string(max_scaled_diff)); + if (max_scaled_diff < param_.outer_aspin_tolerance_) { + report.converged = true; + solution = locally_solved; + return report; + } + + // ----------- Solve global linear system ----------- + // The matrix is kept as is for now. + // TODO: test with proper ASPIN (in what we have now, columns are updated outside the diagonal + // domain-domain interaction blocks in a Gauss-Seidel-like manner), and other variations. + // NOTE: the above should no longer be true, the matrix is ASPIN now (should be). + // +#if 0 + // Compute rhs. + // The modified residual is D * aspin_residual in each domain. + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + for (const auto& domain : domains_) { + auto aspin_res_local = Details::extractVector(aspin_residual, domain.cells); + BVector aspin_res_mod(aspin_res_local.size()); + aspin_res_mod = 0; + if (domain_reports[domain.index].total_newton_iterations > 0) { + (*domain_matrices_[domain.index]).mv(aspin_res_local, aspin_res_mod); + } + Details::setGlobal(aspin_res_mod, domain.cells, ebosResid); + } + // Add J_aspin * aspin_residual to the modified residual, so we solve for + // the difference from the local solution rather than from the very start + // of the iteration. + { + const auto& ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + BVector yy(nc); + ebosJac.mv(aspin_residual, yy); + ebosResid -= yy; + } +#else + // Compute rhs. + { + auto Jcopy = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + Dune::storeMatrixMarket(Jcopy, "aspinjac.mm"); + // Eliminate the domain diagonal blocks from Jcopy. Jcopy is then (J - D) + for (const auto& domain : domains_) { + for (const int cell : domain.cells) { + auto& row = Jcopy[cell]; + for (auto iter = row.begin(); iter != row.end(); ++iter) { + int col = iter.index(); + auto lb = std::find(domain.cells.begin(), domain.cells.end(), col); + if (lb != domain.cells.end()) { + //std::cerr << "row = " << cell << " col = " << col <<'\n'; + *iter = 0.0; + } + } + } + } + Dune::storeMatrixMarket(Jcopy, "aspinjacoffdiag.mm"); + // Set the residual to -Jcopy * aspin_residual, which will then be (D - J)*(x^k - x^{k+1}) + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + ebosResid = 0.0; + Jcopy.mmv(aspin_residual, ebosResid); + } +#endif + + // Check convergence before solving. + { + auto convrep = getConvergence(timer, -iteration, residual_norms); + // report.converged = convrep.converged() && iteration > nonlinear_solver.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.convergence_check_time += perfTimer.stop(); + // residual_norms_history_.push_back(residual_norms); + + // if (report.converged) { + // return report; + // } + + + // Solve the linear system. + linear_solve_setup_time_ = 0.0; + BVector x(nc); + try { + solveJacobianSystem(x); + report.linear_solve_setup_time += linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += linearIterationsLastSolve(); + } + catch (...) { + report.linear_solve_setup_time += linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += linearIterationsLastSolve(); + + failureReport_ += report; + throw; // re-throw up + } + + perfTimer.reset(); + perfTimer.start(); + + // ----------- Update solution ----------- + + // Compute the nonlinear update. + + // apply the Schur compliment of the well model to the reservoir linearized + // equations + // wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), + // ebosSimulator().model().linearizer().residual()); + + + // handling well state update before oscillation treatment is a decision based + // on observation to avoid some big performance degeneration under some circumstances. + // there is no theorectical explanation which way is better for sure. + // wellModel().wellState() = initial_wellstate; + // ebosSimulator_.model().newtonMethod().setIterationIndex(iteration); + // wellModel().beginIteration(); + + // Before we can do postSolve() we must take into account that the current + // state in the wells are from after the local solves, not from the start + // of the iteration. + // Instead of passing x (where state_{i+1} = state_i + x, i.e. x is the + // difference from one global iteration to the next) we should pass + // y, with y = state_{i+1} - state_after_local_solve. + // BVector y(nc); + // for (int c = 0; c < nc; ++c) { + // y[c] = solution[c] - locally_solved[c]; + // } + // wellModel().postSolve(y); + + if (param_.use_update_stabilization_) { + // Stabilize the nonlinear update. + bool isOscillate = false; + bool isStagnate = false; + nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + if (isOscillate) { + current_relaxation_ -= nonlinear_solver.relaxIncrement(); + current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); + if (terminalOutputEnabled()) { + std::string msg = " Oscillating behavior detected: Relaxation set to " + + std::to_string(current_relaxation_); + OpmLog::info(msg); + } + } + nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_); + } + + // Apply the update, with considering model-dependent limitations and + // chopping of the update. This update must be based on the initial solution, + // not the solution after the local solve. + // ebosSimulator().model().solution(0) = initial_solution; + updateSolution(x); + + + + + // Before we can do postSolve() we must take into account that the current + // state in the wells are from after the local solves, not from the start + // of the iteration. + // Instead of passing x (where state_{i+1} = state_i - x, i.e. x is the + // difference from one global iteration to the next) we should pass + // y, with y = state_after_local_solve - state_{i+1}. + // BVector y(nc); + // for (int c = 0; c < nc; ++c) { + // y[c] = locally_solved[c] - solution[c]; + // } + // wellModel().postSolve(y); + wellModel().postSolve(x); + + + + + report.update_time += perfTimer.stop(); + + //wellModel().logPrimaryVars(); + + // Extra debug output. + { + std::ostringstream os; + os << " ** Well rates:"; + for (int w = 0; w < wellModel().wellState().numWells(); ++w) { + os << " | "; + const auto& wr = wellModel().wellState().wellRates(w); + for (const double r : wr) { + os << ' ' << r; + } + } + OpmLog::debug(os.str()); + } + + + return report; + } + + + + + std::pair + solveLocal(const Domain& domain, + const SimulatorTimerInterface& timer, + const int global_iteration, + const bool initial_assembly_required = false) + { + SimulatorReportSingle report; + Dune::Timer solveTimer; + solveTimer.start(); + Dune::Timer detailTimer; + + // TODO: the following assertion will fail, does it indicate a programming error? + // assert(ebosSimulator_.model().newtonMethod().numIterations() == 0); + ebosSimulator_.model().newtonMethod().setIterationIndex(0); + //ebosSimulator_.problem().beginIteration(); + + // 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) { + OpmLog::debug("solveLocal() initial assembly START"); + //wellModel().logPrimaryVars(); + detailTimer.start(); + ebosSimulator_.model().newtonMethod().setIterationIndex(iter); + // TODO: we should have a beginIterationLocal function() + // only handling the well model for now + // ebosSimulator_.problem().beginIteration(); + ebosSimulator_.problem().wellModel().assembleLocal(ebosSimulator_.model().newtonMethod().numIterations(), + ebosSimulator_.timeStepSize(), + domain); + // Assemble reservoir locally. + report += assembleReservoirLocal(domain, iter); + report.assemble_time += detailTimer.stop(); + OpmLog::debug("solveLocal() initial assembly END"); + //wellModel().logPrimaryVars(); + } + detailTimer.reset(); + detailTimer.start(); + std::vector resnorms; + auto convreport = getLocalConvergence(domain, timer, 0, resnorms); + if (convreport.converged()) { + // TODO: set more info, timing etc. + report.converged = true; + // Exit, but not if we are on the first global iteration + // (if so, we require at least one iteration). + //if (global_iteration > 0) { + return { report, convreport }; + //} + } else { + // std::ostringstream os; + // os << "Convergence data for first local iteration:\n" + // << convreport; + // OpmLog::debug(os.str()); + } + //report.convergence_check_time += detailTimer.stop(); + + // 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; + //wellModel().logPrimaryVars(); + + // Local Newton loop. + const int max_iter = param_.max_local_solve_iterations_; + bool converged = false; + 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(); + solveLocalJacobianSystem(domain, x); + wellModel().postSolveLocal(x, domain); + report.linear_solve_time += detailTimer.stop(); + report.linear_solve_setup_time += linear_solve_setup_time_; + report.total_linear_iterations = linearIterationsLastSolve(); + //wellModel().logPrimaryVars(); + + // 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(); + // Note: endIteration does not do anything + // ebosSimulator_.problem().endIteration(); + + // 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 + // ebosSimulator_.problem().beginIteration(); + // Assemble reservoir locally. + ebosSimulator_.problem().wellModel().assembleLocal(ebosSimulator_.model().newtonMethod().numIterations(), + ebosSimulator_.timeStepSize(), + domain); + report += assembleReservoirLocal(domain, iter); + report.assemble_time += detailTimer.stop(); + + + // Check for local convergence. + detailTimer.reset(); + detailTimer.start(); + convreport = getLocalConvergence(domain, timer, iter, resnorms); + // report.convergence_check_time += detailTimer.stop(); + //wellModel().logPrimaryVars(); + // Dune::writeMatrixToMatlab(ebosSimulator_.model().linearizer().jacobian().istlMatrix(), "beforeWellModelLinearize.matlab"); + + // apply the Schur compliment 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; + //wellModel().logPrimaryVars(); + + } while (!convreport.converged() && iter <= max_iter); + + ebosSimulator_.problem().endIteration(); + + if (!convreport.converged()) { + // std::ostringstream os; + // os << "Convergence data for unconverged local solution:\n" + // << convreport; + // OpmLog::debug(os.str()); + } + + 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 }; + } + + + + void printIf(int c, double x, double y, double eps, std::string type) { if (std::abs(x-y) > eps) { std::cout << type << " " <(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.setMatrix(jac); + 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) { - auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); + auto& ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); // set initial guess @@ -619,12 +1578,33 @@ namespace Opm { // account for parallelization properly. since the residual of ECFV // discretizations does not need to be synchronized across processes to be // consistent, this is not relevant for OPM-flow... - ebosSolver.setMatrix(ebosJac); + //ebosSolver.setMatrix(ebosJac); ebosSolver.solve(x); } + /// 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) { @@ -719,7 +1699,8 @@ namespace Opm { /// of the cells associated with a numerical aquifer. std::tuple localConvergenceData(std::vector& R_sum, std::vector& maxCoeff, - std::vector& B_avg) + std::vector& B_avg, + std::vector& maxCoeffCell) { OPM_TIMEBLOCK(localConvergenceData); double pvSumLocal = 0.0; @@ -734,9 +1715,11 @@ namespace Opm { OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + // 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& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); @@ -759,7 +1742,140 @@ namespace Opm { const auto R2 = ebosResid[cell_idx][compIdx]; R_sum[ compIdx ] += R2; - maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue ); + const double Rval = std::abs( R2 ) / pvValue; + if (Rval > maxCoeff[ compIdx ]) { + maxCoeff[ compIdx ] = Rval; + maxCoeffCell[ compIdx ] = cell_idx; + } + } + + if constexpr (has_solvent_) { + B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + const auto R2 = ebosResid[cell_idx][contiSolventEqIdx]; + R_sum[ contiSolventEqIdx ] += R2; + maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_extbo_) { + B_avg[ contiZfracEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiZfracEqIdx]; + R_sum[ contiZfracEqIdx ] += R2; + maxCoeff[ contiZfracEqIdx ] = std::max( maxCoeff[ contiZfracEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_polymer_) { + B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx]; + R_sum[ contiPolymerEqIdx ] += R2; + maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_foam_) { + B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiFoamEqIdx]; + R_sum[ contiFoamEqIdx ] += R2; + maxCoeff[ contiFoamEqIdx ] = std::max( maxCoeff[ contiFoamEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_brine_) { + B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiBrineEqIdx]; + R_sum[ contiBrineEqIdx ] += R2; + maxCoeff[ contiBrineEqIdx ] = std::max( maxCoeff[ contiBrineEqIdx ], std::abs( R2 ) / pvValue ); + } + + if constexpr (has_polymermw_) { + static_assert(has_polymer_); + + B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight + // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations + // TODO: there should be a more general way to determine the scaling-down coefficient + const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.; + R_sum[contiPolymerMWEqIdx] += R2; + maxCoeff[contiPolymerMWEqIdx] = std::max( maxCoeff[contiPolymerMWEqIdx], std::abs( R2 ) / pvValue ); + } + + if constexpr (has_energy_) { + B_avg[ contiEnergyEqIdx ] += 1.0; + const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx]; + R_sum[ contiEnergyEqIdx ] += R2; + maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue ); + } + + } + + 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& R_sum, + std::vector& maxCoeff, + std::vector& B_avg, + std::vector& maxCoeffCell) + { + double pvSumLocal = 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(); + 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); + auto* intQuantsPtr = ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0); + // const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + if (intQuantsPtr == nullptr) { + OpmLog::debug("** no cached intensive quantities for cell " + std::to_string(cell_idx)); + elemCtx.updatePrimaryIntensiveQuantities(0); + intQuantsPtr = ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0); + assert(intQuantsPtr != nullptr); + } + const auto& intQuants = *intQuantsPtr; + + const auto& fs = intQuants.fluidState(); + + const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); + pvSumLocal += pvValue; + + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) + { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + + B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value(); + const auto R2 = ebosResid[cell_idx][compIdx]; + + R_sum[ compIdx ] += R2; + const double Rval = std::abs( R2 ) / pvValue; + if (Rval > maxCoeff[ compIdx ]) { + maxCoeff[ compIdx ] = Rval; + maxCoeffCell[ compIdx ] = cell_idx; + } } if constexpr (has_solvent_) { @@ -842,12 +1958,43 @@ namespace Opm { 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) @@ -870,7 +2017,7 @@ namespace Opm { continue; } elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + // elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); const auto& cellResidual = ebosResid[cell_idx]; @@ -894,6 +2041,119 @@ namespace Opm { return grid_.comm().sum(errorPV); } + ConvergenceReport getLocalReservoirConvergence(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 double pvSum = localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell); + + // compute global sum and max of quantities + // const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal, + // R_sum, maxCoeff, B_avg); + + auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt); + cnvErrorPvFraction /= pvSum; + + 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; + 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 }; + int cell[2] = { -1, maxCoeffCell[compIdx] }; // No cell associated with MB failures. + for (int ii : {0, 1}) { + if (std::isnan(res[ii])) { + report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});//, cell[ii], res[ii]}); + 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});//, cell[ii], res[ii]}); + 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});//, cell[ii], res[ii]}); + 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});//, cell[ii], res[ii]}); + } + } + } + + // 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; + } + + + ConvergenceReport getReservoirConvergence(const double reportTime, const double dt, const int iteration, @@ -906,7 +2166,8 @@ namespace Opm { const int numComp = numEq; Vector R_sum(numComp, 0.0 ); Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); - const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg); + std::vector maxCoeffCell(numComp, -1); + const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell); // compute global sum and max of quantities const auto [ pvSum, numAquiferPvSum ] = @@ -942,24 +2203,25 @@ namespace Opm { CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, CR::ReservoirFailure::Type::Cnv }; double tol[2] = { tol_mb, tol_cnv }; + int cell[2] = { -1, maxCoeffCell[compIdx] }; // No cell associated with MB failures. for (int ii : {0, 1}) { if (std::isnan(res[ii])) { - report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] > maxResidualAllowed()) { - report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] < 0.0) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] > tol[ii]) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); } report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]); } @@ -1001,6 +2263,22 @@ namespace Opm { return report; } + ConvergenceReport getLocalConvergence(const Domain& domain, + const SimulatorTimerInterface& timer, + const int iteration, + std::vector& residual_norms) + { + std::vector B_avg(numEq, 0.0); + auto report = getLocalReservoirConvergence(timer.simulationTimeElapsed(), + timer.currentStepLength(), + iteration, + domain, + B_avg, + residual_norms); + report += wellModel().getLocalWellConvergence(domain, B_avg, false); + return report; + } + /// Compute convergence based on total mass balance (tol_mb) and maximum /// residual mass balance (tol_cnv). /// \param[in] timer simulation timer @@ -1059,6 +2337,10 @@ namespace Opm { const SimulatorReportSingle& failureReport() const { return failureReport_; } + /// return the statistics if the nonlinearIteration() method failed + const SimulatorReportSingle& localAccumulatedReports() const + { return local_reports_accumulated_; } + const std::vector& stepReports() const { return convergence_reports_; @@ -1086,6 +2368,7 @@ namespace Opm { ModelParameters param_; SimulatorReportSingle failureReport_; + SimulatorReportSingle local_reports_accumulated_; // Well Model BlackoilWellModel& well_model_; @@ -1101,6 +2384,9 @@ namespace Opm { std::vector convergence_reports_; ComponentName compNames_{}; + std::vector domains_; + std::vector> domain_matrices_; + std::vector domain_linsolvers_; public: /// return the StandardWells object @@ -1121,6 +2407,19 @@ namespace Opm { } private: + template + struct HasZoltanPartitioning : public std::false_type {}; + + template + struct HasZoltanPartitioning< + GridType, + std::void_t().zoltanPartitionWithoutScatter + (std::declval*>(), + std::declval(), + std::declval(), + std::declval()))> + > : public std::true_type {}; + template bool isNumericalAquiferCell(const Dune::CpGrid& grid, const T& elem) { @@ -1147,6 +2446,46 @@ namespace Opm { double maxResidualAllowed() const { return param_.max_residual_allowed_; } double linear_solve_setup_time_; + std::pair, int> partitionCells() const + { + if constexpr (HasZoltanPartitioning::value) { + return this->partitionCellsZoltan(); + } + else { + return this->partitionCellsFallback(); + } + } + + std::pair, int> partitionCellsZoltan() const + { + const auto wells = this->ebosSimulator_.vanguard().schedule().getWellsatEnd(); + + auto partition_vector = this->grid_.zoltanPartitionWithoutScatter + (&wells, nullptr, this->param_.num_local_domains_, + this->param_.local_domain_partition_imbalance_); + + return this->countDomains(std::move(partition_vector)); + } + + std::pair, int> partitionCellsFallback() const + { + const int num_cells = detail::countLocalInteriorCells(this->grid_); + + return this->countDomains(partitionCells(num_cells)); + } + + std::pair, int> + countDomains(std::vector partition_vector) const + { + auto maxPos = std::max_element(partition_vector.begin(), + partition_vector.end()); + + const auto num_domains = (maxPos == partition_vector.end()) + ? 0 : *maxPos + 1; + + return { std::move(partition_vector), num_domains }; + } + public: std::vector wasSwitched_; }; diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 1da083a476c..62969f9d788 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -168,7 +168,38 @@ template struct UseAverageDensityMsWells { using type = UndefinedProperty; }; - +template +struct NonlinearSolver { + using type = UndefinedProperty; +}; +template +struct LocalSolveApproach { + using type = UndefinedProperty; +}; +template +struct OuterAspinTolerance { + using type = UndefinedProperty; +}; +template +struct MaxLocalSolveIterations { + using type = UndefinedProperty; +}; +template +struct LocalToleranceScalingMb { + using type = UndefinedProperty; +}; +template +struct LocalToleranceScalingCnv { + using type = UndefinedProperty; +}; +template +struct NumLocalDomains { + using type = UndefinedProperty; +}; +template +struct LocalDomainsPartitioningImbalance { + using type = UndefinedProperty; +}; template struct DbhpMaxRel { using type = GetPropType; @@ -316,10 +347,43 @@ template struct UseAverageDensityMsWells { static constexpr bool value = false; }; - - - - +template +struct NonlinearSolver { + static constexpr auto value = "newton"; +}; +template +struct LocalSolveApproach { + static constexpr auto value = "jacobi"; +}; +template +struct OuterAspinTolerance { + using type = GetPropType; + static constexpr type value = 1e-3; +}; +template +struct MaxLocalSolveIterations { + static constexpr int value = 20; +}; +template +struct LocalToleranceScalingMb { + using type = GetPropType; + static constexpr type value = 1.0; +}; +template +struct LocalToleranceScalingCnv { + using type = GetPropType; + static constexpr type value = 0.01; +}; +template +struct NumLocalDomains { + using type = int; + static constexpr auto value = 0; +}; +template +struct LocalDomainsPartitioningImbalance { + using type = GetPropType; + static constexpr auto value = type{1.03}; +}; // if openMP is available, determine the number threads per process automatically. #if _OPENMP @@ -436,7 +500,19 @@ namespace Opm /// Whether to approximate segment densities by averaging over segment and its outlet bool use_average_density_ms_wells_; + // Choose nonlinear solver type: newton, aspin, nldd. + std::string nonlinear_solver_; + std::string local_solve_approach_; + + double outer_aspin_tolerance_; + + int max_local_solve_iterations_; + + double local_tolerance_scaling_mb_; + double local_tolerance_scaling_cnv_; + int num_local_domains_{0}; + double local_domain_partition_imbalance_{1.03}; /// Construct from user parameters or defaults. BlackoilModelParametersEbos() @@ -473,6 +549,14 @@ namespace Opm check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter); max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches); use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells); + nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver); + local_solve_approach_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach); + outer_aspin_tolerance_ = EWOMS_GET_PARAM(TypeTag, double, OuterAspinTolerance); + max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations); + local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb); + local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv); + num_local_domains_ = EWOMS_GET_PARAM(TypeTag, int, NumLocalDomains); + local_domain_partition_imbalance_ = std::max(1.0, EWOMS_GET_PARAM(TypeTag, double, LocalDomainsPartitioningImbalance)); deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); } @@ -512,6 +596,14 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations"); EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, NonlinearSolver, "Choose nonlinear solver. Valid choices are newton, aspin or nldd."); + EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalSolveApproach, "Choose local solve approach. Valid choices are jacobi and gauss-seidel"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, OuterAspinTolerance, "Tolerance for ASPIN residual."); + EWOMS_REGISTER_PARAM(TypeTag, int, MaxLocalSolveIterations, "Max iterations for local solves with ASPIN or NLDD."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingMb, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingCnv, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); + EWOMS_REGISTER_PARAM(TypeTag, int, NumLocalDomains, "Number of local domains for ASPIN/NLDD solves."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalDomainsPartitioningImbalance, "Subdomain partitioning imbalance tolerance. 1.03 is 3 percent imbalance."); } }; } // namespace Opm diff --git a/opm/simulators/flow/FlowMainEbos.hpp b/opm/simulators/flow/FlowMainEbos.hpp index 66c6d2a2057..b24f11892e8 100644 --- a/opm/simulators/flow/FlowMainEbos.hpp +++ b/opm/simulators/flow/FlowMainEbos.hpp @@ -460,6 +460,13 @@ void handleExtraConvergenceOutput(SimulatorReport& report, { SimulatorReport report = simulator_->run(*simtimer_); runSimulatorAfterSim_(report); + { + SimulatorReportSingle rep_local = simulator_->model().localAccumulatedReports(); + std::ostringstream os; + os << "=== Accumulated local solve data ===\n"; + rep_local.reportFullyImplicit(os); + OpmLog::info(os.str()); + } return report.success.exit_status; } diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp index 162f57e2131..efbc724b3c8 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp @@ -362,7 +362,9 @@ class SimulatorFullyImplicitBlackoilEbos // Run a multiple steps of the solver depending on the time step control. solverTimer_->start(); - auto solver = createSolver(wellModel_()); + if (!solver_) { + solver_ = createSolver(wellModel_()); + } ebosSimulator_.startNextEpisode( ebosSimulator_.startTime() @@ -375,8 +377,7 @@ class SimulatorFullyImplicitBlackoilEbos loadStep_ = -1; ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); } - solver->model().beginReportStep(); - + solver_->model().beginReportStep(); bool enableTUNING = EWOMS_GET_PARAM(TypeTag, bool, EnableTuning); // If sub stepping is enabled allow the solver to sub cycle @@ -398,13 +399,13 @@ class SimulatorFullyImplicitBlackoilEbos events.hasEvent(ScheduleEvents::INJECTION_TYPE_CHANGED) || events.hasEvent(ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER) || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE); - auto stepReport = adaptiveTimeStepping_->step(timer, *solver, event, nullptr); + auto stepReport = adaptiveTimeStepping_->step(timer, *solver_, event, nullptr); report_ += stepReport; //Pass simulation report to eclwriter for summary output ebosSimulator_.problem().setSimulationReport(report_); } else { // solve for complete report step - auto stepReport = solver->step(timer); + auto stepReport = solver_->step(timer); report_ += stepReport; if (terminalOutput_) { std::ostringstream ss; @@ -421,7 +422,7 @@ class SimulatorFullyImplicitBlackoilEbos ebosSimulator_.problem().writeOutput(); report_.success.output_write_time += perfTimer.stop(); - solver->model().endReportStep(); + solver_->model().endReportStep(); // take time that was used to solve system for this reportStep solverTimer_->stop(); @@ -433,7 +434,7 @@ class SimulatorFullyImplicitBlackoilEbos // Destructively grab the step convergence reports. The solver // object and the model object contained therein are about to go // out of scope. - this->writeConvergenceOutput(solver->model().getStepReportsDestructively()); + this->writeConvergenceOutput(solver_->model().getStepReportsDestructively()); } // Increment timer, remember well state. @@ -486,6 +487,9 @@ class SimulatorFullyImplicitBlackoilEbos serializer(adaptiveTimeStepping_); } + const Model& model() const + { return solver_->model(); } + protected: std::unique_ptr createSolver(WellModel& wellModel) @@ -702,6 +706,8 @@ class SimulatorFullyImplicitBlackoilEbos ModelParameters modelParam_; SolverParameters solverParam_; + std::unique_ptr solver_; + // Observed objects. PhaseUsage phaseUsage_; // Misc. data diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp new file mode 100644 index 00000000000..db09d823bf4 --- /dev/null +++ b/opm/simulators/flow/SubDomain.hpp @@ -0,0 +1,45 @@ +/* + Copyright 2021 Total SE + + 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_SUBDOMAIN_HEADER_INCLUDED +#define OPM_SUBDOMAIN_HEADER_INCLUDED + +#include + +#include + +namespace Opm +{ + + template + struct SubDomain + { + int index; + std::vector cells; + std::vector interior; + Dune::SubGridPart view; + 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)) + {} + }; + +} // namespace Opm + + +#endif // OPM_SUBDOMAIN_HEADER_INCLUDED diff --git a/opm/simulators/flow/aspinPartition.cpp b/opm/simulators/flow/aspinPartition.cpp new file mode 100644 index 00000000000..6927b721342 --- /dev/null +++ b/opm/simulators/flow/aspinPartition.cpp @@ -0,0 +1,93 @@ +/* + Copyright 2021 Total SE + + 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 + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + +namespace +{ + + // Read from file, containing one number per cell, from [0, ... , num_domains - 1]. + std::pair, int> partitionCellsFromFile([[maybe_unused]] const int num_cells) + { + // TODO: refactor to make more flexible. + // Read file into single vector. + const std::string filename = "partition.txt"; + std::ifstream is(filename); + const std::vector cellpart{std::istream_iterator(is), std::istream_iterator()}; + if (cellpart.size() != size_t(num_cells)) { + auto msg = fmt::format("Partition file contains {} entries, but there are {} cells.", + cellpart.size(), num_cells); + throw std::runtime_error(msg); + } + + // Create and return the output domain vector. + const int num_domains = (*std::max_element(cellpart.begin(), cellpart.end())) + 1; + return { cellpart, num_domains }; + } + + + // Trivially simple partitioner + std::pair, int> partitionCellsSimple(const int num_cells, const int num_domains) + { + // Build the partitions. + std::vector bounds(num_domains + 1); + bounds[0] = 0; + const int dom_sz = num_cells / num_domains; + for (int i = 1; i < num_domains; ++i) { + bounds[i] = dom_sz * i; + } + bounds[num_domains] = num_cells; + std::vector part(num_cells); + for (int i = 0; i < num_domains; ++i) { + std::fill(part.begin() + bounds[i], part.begin() + bounds[i + 1], i); + } + return { part, num_domains }; + } + +} // anonymous namespace + + +std::pair, int> +partitionCells(const int num_cells) +{ + //const std::string method = "simple"; + const std::string method = "file"; + if (method == "simple") { + return partitionCellsSimple(num_cells, 25); + } else if (method == "file") { + return partitionCellsFromFile(num_cells); + } else { + return {}; + } +} + +} // namespace Opm diff --git a/opm/simulators/flow/aspinPartition.hpp b/opm/simulators/flow/aspinPartition.hpp new file mode 100644 index 00000000000..f1781a08696 --- /dev/null +++ b/opm/simulators/flow/aspinPartition.hpp @@ -0,0 +1,34 @@ +/* + Copyright 2021 Total SE + + 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_ASPINPARTITION_HEADER_INCLUDED +#define OPM_ASPINPARTITION_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +std::pair, int> partitionCells(const int num_cells); + +} // namespace Opm + + +#endif // OPM_ASPINPARTITION_HEADER_INCLUDED diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp index 22849fd0f71..3096bafe380 100644 --- a/opm/simulators/linalg/FlexibleSolver.hpp +++ b/opm/simulators/linalg/FlexibleSolver.hpp @@ -101,6 +101,7 @@ class FlexibleSolver : public Dune::InverseOperator preconditioner_; std::shared_ptr scalarproduct_; std::shared_ptr linsolver_; + bool HACK_is_umfpack_ = false; }; } // namespace Dune diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp index 3002c2e25b9..938825c207c 100644 --- a/opm/simulators/linalg/FlexibleSolver_impl.hpp +++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp @@ -73,6 +73,12 @@ namespace Dune FlexibleSolver:: apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) { +#if HAVE_SUITESPARSE_UMFPACK + if (HACK_is_umfpack_) { + using MatrixType = std::remove_const_tgetmat())>>; + linsolver_.reset(new Dune::UMFPack(linearoperator_for_solver_->getmat(), 0, false)); + } +#endif linsolver_->apply(x, rhs, res); } @@ -81,6 +87,12 @@ namespace Dune FlexibleSolver:: apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res) { +#if HAVE_SUITESPARSE_UMFPACK + if (HACK_is_umfpack_) { + using MatrixType = std::remove_const_tgetmat())>>; + linsolver_.reset(new Dune::UMFPack(linearoperator_for_solver_->getmat(), 0, false)); + } +#endif linsolver_->apply(x, rhs, reduction, res); } @@ -185,9 +197,9 @@ namespace Dune verbosity); #if HAVE_SUITESPARSE_UMFPACK } else if (solver_type == "umfpack") { - bool dummy = false; using MatrixType = std::remove_const_tgetmat())>>; - linsolver_ = std::make_shared>(linearoperator_for_solver_->getmat(), verbosity, dummy); + linsolver_ = std::make_shared>(linearoperator_for_solver_->getmat(), verbosity, false); + HACK_is_umfpack_ = true; #endif #if HAVE_CUDA } else if (solver_type == "cubicgstab") { diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 425c41a2c3e..7c7acfa8cfc 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -51,10 +51,6 @@ struct RelaxedLinearSolverReduction { using type = UndefinedProperty; }; -template -struct IluRelaxation { - using type = UndefinedProperty; -}; template struct LinearSolverMaxIter { using type = UndefinedProperty; @@ -67,6 +63,10 @@ struct LinearSolverRestart { // LinearSolverVerbosity defined in opm-models // template +struct IluRelaxation { + using type = UndefinedProperty; +}; +template struct IluFillinLevel { using type = UndefinedProperty; }; @@ -91,23 +91,23 @@ struct LinearSolverIgnoreConvergenceFailure{ using type = UndefinedProperty; }; template -struct PreconditionerAddWellContributions { +struct ScaleLinearSystem { using type = UndefinedProperty; }; template -struct ScaleLinearSystem { +struct LinearSolver { using type = UndefinedProperty; }; template -struct CprReuseSetup { +struct LinearSolverPrintJsonDefinition { using type = UndefinedProperty; }; template -struct CprReuseInterval { +struct CprReuseSetup { using type = UndefinedProperty; }; template -struct LinearSolver { +struct CprReuseInterval { using type = UndefinedProperty; }; template @@ -137,11 +137,6 @@ struct RelaxedLinearSolverReduction { static constexpr type value = 1e-2; }; template -struct IluRelaxation { - using type = GetPropType; - static constexpr type value = 0.9; -}; -template struct LinearSolverMaxIter { static constexpr int value = 200; }; @@ -154,6 +149,11 @@ struct LinearSolverVerbosity { static constexpr int value = 0; }; template +struct IluRelaxation { + using type = GetPropType; + static constexpr type value = 0.9; +}; +template struct IluFillinLevel { static constexpr int value = 0; }; @@ -178,16 +178,16 @@ struct LinearSolverIgnoreConvergenceFailure static constexpr bool value = false; }; template -struct LinearSolverBackend { - using type = ISTLSolverEbos; +struct ScaleLinearSystem { + static constexpr bool value = false; }; template -struct PreconditionerAddWellContributions { - static constexpr bool value = false; +struct LinearSolver { + static constexpr auto value = "ilu0"; }; template -struct ScaleLinearSystem { - static constexpr bool value = false; +struct LinearSolverPrintJsonDefinition { + static constexpr auto value = true; }; template struct CprReuseSetup { @@ -198,10 +198,6 @@ struct CprReuseInterval { static constexpr int value = 30; }; template -struct LinearSolver { - static constexpr auto value = "ilu0"; -}; -template struct AcceleratorMode { static constexpr auto value = "none"; }; @@ -218,6 +214,12 @@ struct OpenclIluParallel { static constexpr bool value = true; // note: false should only be used in debug }; +// Set the backend to be used. +template +struct LinearSolverBackend { + using type = ISTLSolverEbos; +}; + } // namespace Opm::Properties namespace Opm @@ -229,24 +231,24 @@ namespace Opm { double linear_solver_reduction_; double relaxed_linear_solver_reduction_; - double ilu_relaxation_; int linear_solver_maxiter_; int linear_solver_restart_; int linear_solver_verbosity_; + double ilu_relaxation_; int ilu_fillin_level_; MILU_VARIANT ilu_milu_; bool ilu_redblack_; bool ilu_reorder_sphere_; bool newton_use_gmres_; - bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; bool scale_linear_system_; std::string linsolver_; + bool linear_solver_print_json_definition_; + int cpr_reuse_setup_; + int cpr_reuse_interval_; std::string accelerator_mode_; int bda_device_id_; int opencl_platform_id_; - int cpr_reuse_setup_; - int cpr_reuse_interval_; bool opencl_ilu_parallel_; template @@ -255,10 +257,10 @@ namespace Opm // TODO: these parameters have undocumented non-trivial dependencies linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction); relaxed_linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, RelaxedLinearSolverReduction); - ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); linear_solver_maxiter_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter); linear_solver_restart_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverRestart); linear_solver_verbosity_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity); + ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); ilu_fillin_level_ = EWOMS_GET_PARAM(TypeTag, int, IluFillinLevel); ilu_milu_ = convertString2Milu(EWOMS_GET_PARAM(TypeTag, std::string, MiluVariant)); ilu_redblack_ = EWOMS_GET_PARAM(TypeTag, bool, IluRedblack); @@ -266,6 +268,8 @@ namespace Opm newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres); ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure); scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); + linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver); + linear_solver_print_json_definition_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); cpr_reuse_interval_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseInterval); @@ -284,12 +288,12 @@ namespace Opm template static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve should try to achive"); + EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve"); EWOMS_REGISTER_PARAM(TypeTag, double, RelaxedLinearSolverReduction, "The minimum reduction of the residual which the linear solver need to achieve for the solution to be accepted"); - EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIter, "The maximum number of iterations of the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverRestart, "The number of iterations after which GMRES is restarted"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, "The verbosity level of the linear solver (0: off, 2: all)"); + EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, IluFillinLevel, "The fill-in level of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, std::string, MiluVariant, "Specify which variant of the modified-ILU preconditioner ought to be used. Possible variants are: ILU (default, plain ILU), MILU_1 (lump diagonal with dropped row entries), MILU_2 (lump diagonal with the sum of the absolute values of the dropped row entries), MILU_3 (if diagonal is positive add sum of dropped row entrires. Otherwise subtract them), MILU_4 (if diagonal is positive add sum of dropped row entrires. Otherwise do nothing"); EWOMS_REGISTER_PARAM(TypeTag, bool, IluRedblack, "Use red-black partitioning for the ILU preconditioner"); @@ -297,9 +301,10 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge"); EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); + EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition, "Write the JSON definition of the linear solver setup to the DBG wfle."); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval"); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter."); - EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); EWOMS_REGISTER_PARAM(TypeTag, std::string, AcceleratorMode, "Choose a linear solver, usage: '--accelerator-mode=[none|cusparse|opencl|amgcl|rocalution]'"); EWOMS_REGISTER_PARAM(TypeTag, int, BdaDeviceId, "Choose device ID for cusparseSolver or openclSolver, use 'nvidia-smi' or 'clinfo' to determine valid IDs"); EWOMS_REGISTER_PARAM(TypeTag, int, OpenclPlatformId, "Choose platform ID for openclSolver, use 'clinfo' to determine valid platform IDs"); @@ -311,19 +316,23 @@ namespace Opm // set default values void reset() { - newton_use_gmres_ = false; - linear_solver_reduction_ = 1e-2; relaxed_linear_solver_reduction_ = 1e-2; - linear_solver_maxiter_ = 150; - linear_solver_restart_ = 40; - linear_solver_verbosity_ = 0; - require_full_sparsity_pattern_ = false; - ignoreConvergenceFailure_ = false; - ilu_fillin_level_ = 0; + linear_solver_reduction_ = 1e-2; + linear_solver_maxiter_ = 200; + linear_solver_restart_ = 40; + linear_solver_verbosity_ = 0; ilu_relaxation_ = 0.9; + ilu_fillin_level_ = 0; ilu_milu_ = MILU_VARIANT::ILU; ilu_redblack_ = false; - ilu_reorder_sphere_ = true; + ilu_reorder_sphere_ = false; + newton_use_gmres_ = false; + ignoreConvergenceFailure_ = false; + scale_linear_system_ = false; + linsolver_ = "ilu0"; + linear_solver_print_json_definition_ = true; + cpr_reuse_setup_ = 4; + cpr_reuse_interval_ = 30; accelerator_mode_ = "none"; bda_device_id_ = 0; opencl_platform_id_ = 0; diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 068a3ce93e4..f4d9194ea97 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -231,25 +231,45 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } /// Construct a system solver. - /// \param[in] parallelInformation In the case of a parallel run - /// with dune-istl the information about the parallelization. + /// \param[in] simulator The opm-models simulator object + /// \param[in] parameters Explicit parameters for solver setup, do not + /// read them from command line parameters. + ISTLSolverEbos(const Simulator& simulator, const FlowLinearSolverParameters& parameters) + : simulator_(simulator), + iterations_( 0 ), + calls_( 0 ), + converged_(false), + matrix_(nullptr), + parameters_(parameters) + { + initialize(); + } + + /// Construct a system solver. + /// \param[in] simulator The opm-models simulator object explicit ISTLSolverEbos(const Simulator& simulator) : simulator_(simulator), iterations_( 0 ), calls_( 0 ), converged_(false), - matrix_() + matrix_(nullptr) { - OPM_TIMEBLOCK(IstlSolverEbos); - const bool on_io_rank = (simulator.gridView().comm().rank() == 0); -#if HAVE_MPI - comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); -#endif parameters_.template init(simulator_.vanguard().eclState().getSimulationConfig().useCPR()); + initialize(); + } + + void initialize() + { + OPM_TIMEBLOCK(IstlSolverEbos); prm_ = setupPropertyTree(parameters_, EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter), EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)); + const bool on_io_rank = (simulator_.gridView().comm().rank() == 0); +#if HAVE_MPI + comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); +#endif + #if COMPILE_BDA_BRIDGE { std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode); @@ -300,7 +320,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true); // Print parameters to PRT/DBG logs. - if (on_io_rank) { + if (on_io_rank && parameters_.linear_solver_print_json_definition_) { std::ostringstream os; os << "Property tree for linear solver:\n"; prm_.write_json(os, true); @@ -314,12 +334,17 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } void prepare(const SparseMatrixAdapter& M, Vector& b) + { + prepare(M.istlMatrix(), b); + } + + void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverEbosPrepare); - static bool firstcall = true; + const bool firstcall = (matrix_ == nullptr); #if HAVE_MPI if (firstcall) { - const size_t size = M.istlMatrix().N(); + const size_t size = M.N(); detail::copyParValues(parallelInformation_, size, *comm_); } #endif @@ -329,7 +354,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, // ebos will not change the matrix object. Hence simply store a pointer // to the original one with a deleter that does nothing. // Outch! We need to be able to scale the linear system! Hence const_cast - matrix_ = const_cast(&M.istlMatrix()); + matrix_ = const_cast(&M); useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver) @@ -343,7 +368,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #endif } else { // Pointers should not change - if ( &(M.istlMatrix()) != matrix_ ) { + if ( &M != matrix_ ) { OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"); } @@ -360,8 +385,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #else prepareFlexibleSolver(); #endif - - firstcall = false; } diff --git a/opm/simulators/linalg/extractMatrix.hpp b/opm/simulators/linalg/extractMatrix.hpp new file mode 100644 index 00000000000..b659f2e28b3 --- /dev/null +++ b/opm/simulators/linalg/extractMatrix.hpp @@ -0,0 +1,151 @@ +/* + Copyright 2021 SINTEF Digital, Mathematics and Cybernetics. + + 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_EXTRACT_MATRIX_HEADER_INCLUDED +#define OPM_EXTRACT_MATRIX_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +namespace Details +{ + + template + void copySubMatrix(const Matrix& A, const std::vector& indices, Matrix& B) + { + // Copy elements so that B_{i, j} = A_{indices[i], indices[j]}. + for (auto row = B.begin(); row != B.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + *col = A[indices[row.index()]][indices[col.index()]]; + } + } + } + + + template + Matrix extractMatrix(const Matrix& m, const std::vector& indices) + { + assert(std::is_sorted(indices.begin(), indices.end())); + + // Set up reverse index map. + const size_t n = indices.size(); + size_t newrow = 0; + enum { NotIncluded = -1 }; + std::vector index_map(m.N(), NotIncluded); + for (auto row = m.begin(); row != m.end(); ++row) { + if (row.index() == indices[newrow]) { + index_map[row.index()] = newrow; + ++newrow; + if (newrow == n) { + break; + } + } + } + assert(newrow == n); + + // Count nonzeroes. + size_t nnz = 0; + for (auto row = m.begin(); row != m.end(); ++row) { + if (index_map[row.index()] != NotIncluded) { + for (auto col = row->begin(); col != row->end(); ++col) { + if (index_map[col.index()] != NotIncluded) { + ++nnz; + } + } + } + } + + // Create the matrix structure. + Matrix res(n, n, nnz, Matrix::row_wise); + auto from_row = m.begin(); + for (auto row = res.createbegin(); row != res.createend(); ++row) { + // Move from_row to point to the row to be extracted. + while (from_row.index() < indices[row.index()]) { + ++from_row; + } + assert(from_row.index() == indices[row.index()]); + // Insert nonzeros for row. + for (auto from_col = from_row->begin(); from_col != from_row->end(); ++from_col) { + const int new_col = index_map[from_col.index()]; + if (new_col != NotIncluded) { + row.insert(new_col); + } + } + } + + copySubMatrix(m, indices, res); + return res; + } + + + template + Vector extractVector(const Vector& x, const std::vector& indices) + { + Vector res(indices.size()); + for (size_t ii = 0; ii < indices.size(); ++ii) { + res[ii] = x[indices[ii]]; + } + return res; + } + + + template + void setGlobal(const Vector& x, const std::vector& indices, Vector& global_x) + { + for (size_t ii = 0; ii < indices.size(); ++ii) { + global_x[indices[ii]] = x[ii]; + } + } + + + + template + bool matrixEqual(const Matrix& m1, const Matrix& m2) + { + // Compare size and nonzeroes. + if (m1.N() != m2.N()) return false; + if (m1.M() != m2.M()) return false; + if (m1.nonzeroes() != m2.nonzeroes()) return false; + + auto row1 = m1.begin(); + auto row2 = m2.begin(); + for (; row1 != m1.end(); ++row1, ++row2) { + if (row2 == m2.end()) return false; + if (row1.index() != row2.index()) return false; + auto col1 = row1->begin(); + auto col2 = row2->begin(); + for (; col1 != row1->end(); ++col1, ++col2) { + if (col2 == row2->end()) return false; + if (col1.index() != col2.index()) return false; + if (*col1 != *col2) return false; + } + } + return true; + } + + +} // namespace Details + + +} // namespace Opm + +#endif // OPM_EXTRACT_MATRIX_HEADER_INCLUDED diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 0abbaad9d2a..6b731040e19 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -45,6 +45,7 @@ #include #include +#include #include #include #include @@ -137,6 +138,8 @@ namespace Opm { using AverageRegionalPressureType = RegionAverageCalculator:: AverageRegionalPressure >; + using Domain = SubDomain; + BlackoilWellModel(Simulator& ebosSimulator); void init(); @@ -155,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 postSolveLocal(GlobalEqVector& deltaX, const Domain& domain) + { + recoverWellSolutionAndUpdateWellStateLocal(deltaX, domain); + } + ///////////// // ///////////// @@ -270,6 +279,11 @@ 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 getLocalWellConvergence(const Domain& domain, + const std::vector& B_avg, + const bool checkGroupConvergence = false) const; + const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const; @@ -292,6 +306,7 @@ namespace Opm { // at the beginning of each time step (Not report step) void prepareTimeStep(DeferredLogger& deferred_logger); void initPrimaryVariablesEvaluation() const; + void initPrimaryVariablesEvaluationLocal(const Domain& domain) const; std::pair updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false); @@ -328,6 +343,20 @@ namespace Opm { int numLocalNonshutWells() const; + // prototype for assemble function for ASPIN solveLocal() + // will try to merge back to assemble() when done prototyping + void assembleLocal(const int iterationIdx, + const double dt, + const Domain& domain); + void updateWellControlsLocal(DeferredLogger& deferred_logger, const bool checkGroupControls, + const Domain& domain); + + void logPrimaryVars() const; + std::vector getPrimaryVarsDomain(const Domain& domain) const; + void setPrimaryVarsDomain(const Domain& domain, const std::vector& vars); + + void setupDomains(const std::vector& domains); + protected: Simulator& ebosSimulator_; @@ -377,6 +406,9 @@ namespace Opm { std::vector B_avg_{}; + // Keep track of the domain of each well, if using subdomains. + std::map well_domain_; + const Grid& grid() const { return ebosSimulator_.vanguard().grid(); } @@ -404,6 +436,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); @@ -414,6 +447,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 recoverWellSolutionAndUpdateWellStateLocal(const BVector& x, const Domain& domain); + // setting the well_solutions_ based on well_state. void updatePrimaryVariables(DeferredLogger& deferred_logger); @@ -438,6 +475,7 @@ namespace Opm { int reportStepIndex() const; void assembleWellEq(const double dt, DeferredLogger& deferred_logger); + void assembleWellEqLocal(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 aa47b599e59..aa5012ce7c2 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -38,6 +38,7 @@ #endif #include +#include #include #include @@ -169,6 +170,36 @@ namespace Opm { } + template + void + BlackoilWellModel:: + linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res) + { + if (!param_.matrix_add_well_contributions_) { + // 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_) { + if (well_domain_.at(well->name()) == domain.index) { + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + } + return; + } + + for (const auto& well: well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->addWellContributions(jacobian); + // applying the well residual to reservoir residuals + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + } + } + + + + template void BlackoilWellModel:: @@ -1026,6 +1057,66 @@ namespace Opm { + template + void + BlackoilWellModel:: + assembleLocal(const int iterationIdx, + const double dt, + const Domain& domain) + { + + DeferredLogger local_deferredLogger; + // TODO: not handling gas lift for now + /* if (this->glift_debug) { + const std::string msg = fmt::format( + "assemble() : iteration {}" , iterationIdx); + gliftDebug(msg, local_deferredLogger); + } */ + last_report_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + perfTimer.start(); + + /* if ( ! wellsActive() ) { + return; + } */ + + + // TODO: this function remains to be optimized, because we do not need to handle all the wells + updatePerforationIntensiveQuantities(); + + auto exc_type = ExceptionType::NONE; + std::string exc_msg; + try { + // if (iterationIdx == 0) { + // calculateExplicitQuantities(local_deferredLogger); + // prepareTimeStep(local_deferredLogger); + // } + updateWellControlsLocal(local_deferredLogger, /* check group controls */ false, domain); + + // Set the well primary variables based on the value of well solutions + initPrimaryVariablesEvaluationLocal(domain); + + // maybeDoGasLiftOptimize(local_deferredLogger); + assembleWellEqLocal(dt, domain, local_deferredLogger); + } catch (const std::runtime_error& e) { + exc_type = ExceptionType::RUNTIME_ERROR; + exc_msg = e.what(); + } catch (const std::invalid_argument& e) { + exc_type = ExceptionType::INVALID_ARGUMENT; + exc_msg = e.what(); + } catch (const std::logic_error& e) { + exc_type = ExceptionType::LOGIC_ERROR; + exc_msg = e.what(); + } catch (const std::exception& e) { + exc_type = ExceptionType::DEFAULT; + exc_msg = e.what(); + } + // @TODO errors here must be caught higher up, as this method is not called in parallel. + logAndCheckForExceptionsAndThrow(local_deferredLogger, exc_type, "assemble() failed: " + exc_msg, terminal_output_, Parallel::Communication()); + last_report_.converged = true; + last_report_.assemble_time_well += perfTimer.stop(); + } + template bool @@ -1237,6 +1328,17 @@ namespace Opm { well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger); } } + template + void + BlackoilWellModel:: + assembleWellEqLocal(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 @@ -1466,7 +1568,6 @@ namespace Opm { BlackoilWellModel:: recoverWellSolutionAndUpdateWellState(const BVector& x) { - DeferredLogger local_deferredLogger; OPM_BEGIN_PARALLEL_TRY_CATCH(); { @@ -1474,12 +1575,34 @@ namespace Opm { 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:: + recoverWellSolutionAndUpdateWellStateLocal(const BVector& x, const Domain& domain) + { + DeferredLogger local_deferredLogger; + OPM_BEGIN_PARALLEL_TRY_CATCH(); + { + 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); + } + } + } + OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, + "recoverWellSolutionAndUpdateWellState() failed: ", + terminal_output_, ebosSimulator_.vanguard().grid().comm()); } @@ -1495,6 +1618,81 @@ namespace Opm { } } + template + void + BlackoilWellModel:: + initPrimaryVariablesEvaluationLocal(const Domain& domain) const + { + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->initPrimaryVariablesEvaluation(); + } + } + } + + + + + + + template + ConvergenceReport + BlackoilWellModel:: + getLocalWellConvergence(const Domain& domain, + const std::vector& B_avg, + const bool checkWellGroupControls) const + { + const auto& summary_state = ebosSimulator_.vanguard().summaryState(); + + const auto relax_tolerance = ebosSimulator_.model().newtonMethod().numIterations() + > param_.strict_outer_iter_wells_; + + Opm::DeferredLogger local_deferredLogger; + + // Get global (from all processes) convergence report. + ConvergenceReport local_report; + + for (const auto& well : well_container_) { + if ((well_domain_.at(well->name()) == domain.index) && + well->isOperableAndSolvable()) + { + local_report += well->getWellConvergence(summary_state, this->wellState(), B_avg, + local_deferredLogger, relax_tolerance); + } + } + + // 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. + Opm::DeferredLogger global_deferredLogger = local_deferredLogger; + ConvergenceReport report = local_report; + + // the well_group_control_changed info is already communicated + // TODO: Check what the implications are here for local solves. + if (checkWellGroupControls) { + report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed); + } + + if (terminal_output_) { + global_deferredLogger.logMessages(); + } + + + // Log debug messages for NaN or too large residuals. + 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 +1839,55 @@ namespace Opm { return { changed_well_group, more_network_update }; } + template + void + BlackoilWellModel:: + updateWellControlsLocal(DeferredLogger& deferred_logger, const bool checkGroupControls, const Domain& domain) + { + // Even if there are no wells active locally, we cannot + // return as the DeferredLogger uses global communication. + // For no well active globally we simply return. + if( !wellsActive() ) return ; + +/* updateAndCommunicateGroupData(); + + updateNetworkPressures(); */ + + std::set switched_wells; + std::set switched_groups; + +/* if (checkGroupControls) { + // Check group individual constraints. + updateGroupIndividualControls(deferred_logger, switched_groups); + + // Check group's constraints from higher levels. + updateGroupHigherControls(deferred_logger, switched_groups); + + updateAndCommunicateGroupData(); + + // Check wells' group constraints and communicate. + for (const auto& well : well_container_) { + const auto mode = WellInterface::IndividualOrGroup::Group; + const bool changed = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); + if (changed) { + switched_wells.insert(well->name()); + } + } + updateAndCommunicateGroupData(); + } */ + + // Check individual well constraints and communicate. + for (const auto& well : well_container_) { + if (switched_wells.count(well->name()) || well_domain_.at(well->name()) != domain.index) { + continue; + } + const auto mode = WellInterface::IndividualOrGroup::Individual; + well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); + } + // updateAndCommunicateGroupData(); + } + + template void @@ -1949,6 +2196,7 @@ namespace Opm { BlackoilWellModel:: updatePerforationIntensiveQuantities() { +#if 0 ElementContext elemCtx(ebosSimulator_); const auto& gridView = ebosSimulator_.gridView(); @@ -1963,6 +2211,7 @@ namespace Opm { elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); } OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updatePerforationIntensiveQuantities() failed: ", ebosSimulator_.vanguard().grid().comm()); +#endif } @@ -2094,6 +2343,59 @@ namespace Opm { } + template + void + BlackoilWellModel:: + logPrimaryVars() const + { + std::ostringstream os; + for (const auto& w : well_container_) { + os << w->name() << ":"; + auto pv = w->getPrimaryVars(); + for (const double v : pv) { + os << ' ' << v; + } + os << '\n'; + } + OpmLog::debug(os.str()); + } + + + + template + std::vector + BlackoilWellModel:: + getPrimaryVarsDomain(const Domain& domain) const + { + std::vector ret; + for (const auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + const auto& pv = well->getPrimaryVars(); + for (const double v : pv) { + ret.push_back(v); + } + } + } + return ret; + } + + + + template + void + BlackoilWellModel:: + setPrimaryVarsDomain(const Domain& domain, const std::vector& vars) + { + int offset = 0; + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + int num_pri_vars = well->setPrimaryVars(vars.begin() + offset); + offset += num_pri_vars; + } + } + } + + template void @@ -2118,7 +2420,38 @@ namespace Opm { } - - + template + void + BlackoilWellModel:: + setupDomains(const std::vector& domains) + { + for (const auto& wellPtr : this->well_container_) { + const int first_well_cell = wellPtr->cells()[0]; + for (const auto& domain : domains) { + const bool found = std::binary_search(domain.cells.begin(), domain.cells.end(), first_well_cell); + if (found) { + // Assuming that if the first well cell is found in a domain, + // then all of that well's cells are in that same domain. + // TODO verify assumption. + 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 (!std::binary_search(domain.cells.begin(), domain.cells.end(), well_cell)) { + // TODO: ensure proper exit in parallel. + OPM_THROW(std::runtime_error, "Well found on multiple domains."); + } + } + } + } + } + // TODO: ensure output on rank 0 only. + std::ostringstream os; + os << "Well name Domain\n"; + for (const auto& [wname, domain] : well_domain_) { + os << wname << std::setw(21 - wname.size()) << domain << '\n'; + } + OpmLog::debug(os.str()); + } } // namespace Opm diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index d2febff4af1..4c48699a497 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -254,6 +254,25 @@ namespace Opm DeferredLogger& deferred_logger) const; + std::vector getPrimaryVars() const override + { + const int num_pri_vars = this->primary_variables_.numWellEq(); + std::vector retval(num_pri_vars); + for (int ii = 0; ii < num_pri_vars; ++ii) { + retval[ii] = this->primary_variables_.value(ii); + } + return retval; + } + + int setPrimaryVars(std::vector::const_iterator it) override + { + const int num_pri_vars = this->primary_variables_.numWellEq(); + for (int ii = 0; ii < num_pri_vars; ++ii) { + this->primary_variables_.setValue(ii, it[ii]); + } + return num_pri_vars; + } + protected: bool regularize_; diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 31951536ede..10037f6771a 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -135,7 +135,7 @@ class StandardWellPrimaryVariables { //! \brief Returns scaled rate for a component. EvalWell getQs(const int compIdx) const; - //! \brief Returns a const ref to an evaluation. + //! \brief Returns a value. Scalar value(const int idx) const { return value_[idx]; } @@ -143,6 +143,10 @@ class StandardWellPrimaryVariables { const EvalWell& eval(const int idx) const { return evaluation_[idx]; } + //! \brief Set a value. Note that this does not also set the corresponding evaluation. + void setValue(const int idx, const Scalar val) + { value_[idx] = val; } + private: //! \brief Calculate a relaxation factor for producers. //! \details To avoid overshoot of the fractions which might result in negative rates. diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index d3a2fa336ad..cbb49126fed 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -328,6 +328,16 @@ class WellInterface : public WellInterfaceIndices getPrimaryVars() const + { + return {}; + } + + virtual int setPrimaryVars(std::vector::const_iterator) + { + return 0; + } + protected: // simulation parameters const ModelParameters& param_; diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 754372b8c04..88c6e5fc135 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -397,6 +397,11 @@ bool WellInterfaceGeneric::isVFPActive(DeferredLogger& deferred_logger) const } } +const Well& WellInterfaceGeneric::scheduleWell() const +{ + return well_ecl_; +} + bool WellInterfaceGeneric::isOperableAndSolvable() const { return operability_status_.isOperableAndSolvable(); diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 674e50fff27..852b0956de8 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -197,6 +197,9 @@ class WellInterfaceGeneric { bool stopppedOrZeroRateTarget(const SummaryState& summary_state, const WellState& well_state) const; + + const Well& scheduleWell() const; + protected: bool getAllowCrossFlow() const; diff --git a/tests/partition.txt b/tests/partition.txt new file mode 100644 index 00000000000..86d2c9c5f94 --- /dev/null +++ b/tests/partition.txt @@ -0,0 +1,10 @@ +0 +0 +1 +1 +2 +2 +1 +1 +0 +0 diff --git a/tests/test_aspinPartition.cpp b/tests/test_aspinPartition.cpp new file mode 100644 index 00000000000..de542d81c52 --- /dev/null +++ b/tests/test_aspinPartition.cpp @@ -0,0 +1,35 @@ +/* + Copyright 2021 SINTEF Digital, Mathematics and Cybernetics. + + 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 + +#define BOOST_TEST_MODULE OPM_test_aspinPartition +#include + +#include + + +BOOST_AUTO_TEST_CASE(fileBased) +{ + auto [part, num_part] = Opm::partitionCells(10); + BOOST_CHECK_EQUAL(num_part, 3); + std::vector expected = { 0, 0, 1, 1, 2, 2, 1, 1, 0, 0 }; + BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end()); +} +