82 using BVector =
typename BlackoilModelEbos<TypeTag>::BVector;
85 using Mat =
typename BlackoilModelEbos<TypeTag>::Mat;
87 static constexpr int numEq = Indices::numEq;
96 const auto& grid = model_.ebosSimulator().vanguard().grid();
97 const auto& schedule = model_.ebosSimulator().vanguard().schedule();
102 schedule.getWellsatEnd(),
103 model_.param().local_domain_partition_method_,
104 model_.param().num_local_domains_,
105 model_.param().local_domain_partition_imbalance_);
125 const auto&
gridView = grid.leafGridView();
129 for (
auto it =
beg; it != end; ++it, ++cell) {
138 for (
int index = 0; index <
num_domains; ++index) {
144 Dune::SubGridPart<Grid> view{grid, std::move(
seeds[index])};
146 this->domains_.emplace_back(index,
156 for (
int index = 0; index <
num_domains; ++index) {
160 const auto& eclState = model_.ebosSimulator().vanguard().eclState();
165 if (domains_[index].cells.size() < 200) {
171 loc_param.linear_solver_print_json_definition_ =
false;
183 model_.wellModel().setupDomains(domains_);
187 template <
class NonlinearSolverType>
198 if (report.converged) {
204 auto&
solution = model_.ebosSimulator().model().solution(0);
212 std::vector<SimulatorReportSingle>
domain_reports(domains_.size());
216 switch (model_.param().local_solve_approach_) {
217 case DomainSolveApproach::Jacobi:
222 case DomainSolveApproach::GaussSeidel:
232 OpmLog::debug(
"Convergence failure in domain " + std::to_string(
domain.index));
247 std::ostringstream
os;
248 os << fmt::format(
"Local solves finished. Converged for {}/{} domains.\n",
250 rep.reportFullyImplicit(
os,
nullptr);
251 OpmLog::debug(
os.str());
252 local_reports_accumulated_ +=
rep;
255 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
257 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(0);
268 const auto& comm = model_.ebosSimulator().vanguard().grid().comm();
269 if (comm.size() > 1) {
270 const auto*
ccomm = model_.ebosSimulator().model().newtonMethod().linearSolver().comm();
278 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
282 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
287 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(0);
299 report.converged =
true;
307 return local_reports_accumulated_;
312 std::pair<SimulatorReportSingle, ConvergenceReport>
313 solveDomain(
const Domain&
domain,
318 auto& ebosSimulator = model_.ebosSimulator();
325 ebosSimulator.model().newtonMethod().setIterationIndex(0);
334 ebosSimulator.model().newtonMethod().setIterationIndex(
iter);
337 ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(),
338 ebosSimulator.timeStepSize(),
341 report += this->assembleReservoirDomain(
domain);
346 std::vector<double> resnorms;
347 auto convreport = this->getDomainConvergence(domain, timer, 0, resnorms);
348 if (convreport.converged()) {
350 report.converged =
true;
351 return { report, convreport };
358 model_.wellModel().linearizeDomain(domain,
359 ebosSimulator.model().linearizer().jacobian(),
360 ebosSimulator.model().linearizer().residual());
361 const double tt1 = detailTimer.stop();
362 report.assemble_time += tt1;
363 report.assemble_time_well += tt1;
366 const int max_iter = model_.param().max_local_solve_iterations_;
367 const auto& grid = ebosSimulator.vanguard().grid();
371 const int nc = grid.size(0);
375 this->solveJacobianSystemDomain(domain, x);
376 model_.wellModel().postSolveDomain(x, domain);
377 report.linear_solve_time += detailTimer.stop();
378 report.linear_solve_setup_time += model_.linearSolveSetupTime();
379 report.total_linear_iterations = model_.linearIterationsLastSolve();
384 this->updateDomainSolution(domain, x);
385 report.update_time += detailTimer.stop();
391 ebosSimulator.model().newtonMethod().setIterationIndex(iter);
395 ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(),
396 ebosSimulator.timeStepSize(),
398 report += this->assembleReservoirDomain(domain);
399 report.assemble_time += detailTimer.stop();
404 convreport = this->getDomainConvergence(domain, timer, iter, resnorms);
410 model_.wellModel().linearizeDomain(domain,
411 ebosSimulator.model().linearizer().jacobian(),
412 ebosSimulator.model().linearizer().residual());
413 const double tt2 = detailTimer.stop();
414 report.assemble_time += tt2;
415 report.assemble_time_well += tt2;
416 }
while (!convreport.converged() && iter <= max_iter);
418 ebosSimulator.problem().endIteration();
420 report.converged = convreport.converged();
421 report.total_newton_iterations = iter;
422 report.total_linearizations = iter;
423 report.total_time = solveTimer.stop();
425 return { report, convreport };
429 SimulatorReportSingle assembleReservoirDomain(
const Domain& domain)
432 model_.ebosSimulator().model().linearizer().linearizeDomain(domain);
433 return model_.wellModel().lastReport();
437 void solveJacobianSystemDomain(
const Domain& domain, BVector& global_x)
439 const auto& ebosSimulator = model_.ebosSimulator();
441 Dune::Timer perfTimer;
444 const Mat& main_matrix = ebosSimulator.model().linearizer().jacobian().istlMatrix();
445 if (domain_matrices_[domain.index]) {
446 Details::copySubMatrix(main_matrix, domain.cells, *domain_matrices_[domain.index]);
448 domain_matrices_[domain.index] = std::make_unique<Mat>(Details::extractMatrix(main_matrix, domain.cells));
450 auto& jac = *domain_matrices_[domain.index];
451 auto res = Details::extractVector(ebosSimulator.model().linearizer().residual(),
459 auto& linsolver = domain_linsolvers_[domain.index];
461 linsolver.prepare(jac, res);
462 model_.linearSolveSetupTime() = perfTimer.stop();
463 linsolver.setResidual(res);
466 Details::setGlobal(x, domain.cells, global_x);
470 void updateDomainSolution(
const Domain& domain,
const BVector& dx)
472 auto& ebosSimulator = model_.ebosSimulator();
473 auto& ebosNewtonMethod = ebosSimulator.model().newtonMethod();
474 SolutionVector& solution = ebosSimulator.model().solution(0);
476 ebosNewtonMethod.update_(solution,
485 ebosSimulator.model().invalidateAndUpdateIntensiveQuantities(0, domain);
489 std::pair<double, double> localDomainConvergenceData(
const Domain& domain,
490 std::vector<Scalar>& R_sum,
491 std::vector<Scalar>& maxCoeff,
492 std::vector<Scalar>& B_avg,
493 std::vector<int>& maxCoeffCell)
495 const auto& ebosSimulator = model_.ebosSimulator();
497 double pvSumLocal = 0.0;
498 double numAquiferPvSumLocal = 0.0;
499 const auto& ebosModel = ebosSimulator.model();
500 const auto& ebosProblem = ebosSimulator.problem();
502 const auto& ebosResid = ebosSimulator.model().linearizer().residual();
504 ElementContext elemCtx(ebosSimulator);
505 const auto& gridView = domain.view;
506 const auto& elemEndIt = gridView.template end<0>();
507 IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
509 for (
auto elemIt = gridView.template begin</*codim=*/0>();
513 if (elemIt->partitionType() != Dune::InteriorEntity) {
516 const auto& elem = *elemIt;
517 elemCtx.updatePrimaryStencil(elem);
518 elemCtx.updatePrimaryIntensiveQuantities(0);
520 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
521 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
522 const auto& fs = intQuants.fluidState();
524 const auto pvValue = ebosProblem.referencePorosity(cell_idx, 0) *
525 ebosModel.dofTotalVolume(cell_idx);
526 pvSumLocal += pvValue;
528 if (isNumericalAquiferCell(elem))
530 numAquiferPvSumLocal += pvValue;
533 model_.getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue,
534 B_avg, R_sum, maxCoeff, maxCoeffCell);
538 const int bSize = B_avg.size();
539 for (
int i = 0; i<bSize; ++i )
541 B_avg[ i ] /= Scalar(domain.cells.size());
544 return {pvSumLocal, numAquiferPvSumLocal};
547 ConvergenceReport getDomainReservoirConvergence(
const double reportTime,
550 const Domain& domain,
551 std::vector<Scalar>& B_avg,
552 std::vector<Scalar>& residual_norms)
554 using Vector = std::vector<Scalar>;
556 const int numComp = numEq;
557 Vector R_sum(numComp, 0.0 );
558 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest() );
559 std::vector<int> maxCoeffCell(numComp, -1);
560 const auto [ pvSum, numAquiferPvSum]
561 = this->localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell);
563 auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt);
564 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
566 const double tol_mb = model_.param().local_tolerance_scaling_mb_ *
567 model_.param().tolerance_mb_;
571 const bool use_relaxed = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
572 iteration >= model_.param().min_strict_cnv_iter_;
575 const double tol_cnv = model_.param().local_tolerance_scaling_cnv_
576 * (use_relaxed ? model_.param().tolerance_cnv_relaxed_
577 : model_.param().tolerance_cnv_);
580 std::vector<Scalar> CNV(numComp);
581 std::vector<Scalar> mass_balance_residual(numComp);
582 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
584 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
585 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
586 residual_norms.push_back(CNV[compIdx]);
590 ConvergenceReport report{reportTime};
591 using CR = ConvergenceReport;
592 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
593 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
594 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
595 CR::ReservoirFailure::Type::Cnv };
596 double tol[2] = { tol_mb, tol_cnv };
597 for (
int ii : {0, 1}) {
598 if (std::isnan(res[ii])) {
599 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
600 if (model_.terminalOutputEnabled()) {
601 OpmLog::debug(
"NaN residual for " + model_.compNames().name(compIdx) +
" equation.");
603 }
else if (res[ii] > model_.param().max_residual_allowed_) {
604 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
605 if (model_.terminalOutputEnabled()) {
606 OpmLog::debug(
"Too large residual for " + model_.compNames().name(compIdx) +
" equation.");
608 }
else if (res[ii] < 0.0) {
609 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
610 if (model_.terminalOutputEnabled()) {
611 OpmLog::debug(
"Negative residual for " + model_.compNames().name(compIdx) +
" equation.");
613 }
else if (res[ii] > tol[ii]) {
614 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
620 if (model_.terminalOutputEnabled())
623 if (iteration == 0) {
624 std::string msg = fmt::format(
"Domain {}, size {}, containing cell {}\n| Iter",
625 domain.index, domain.cells.size(), domain.cells[0]);
626 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
628 msg += model_.compNames().name(compIdx)[0];
631 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
633 msg += model_.compNames().name(compIdx)[0];
638 std::ostringstream ss;
640 const std::streamsize oprec = ss.precision(3);
641 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
642 ss << std::setw(4) << iteration;
643 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
644 ss << std::setw(11) << mass_balance_residual[compIdx];
646 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
647 ss << std::setw(11) << CNV[compIdx];
651 OpmLog::debug(ss.str());
657 ConvergenceReport getDomainConvergence(
const Domain& domain,
658 const SimulatorTimerInterface& timer,
660 std::vector<double>& residual_norms)
662 std::vector<Scalar> B_avg(numEq, 0.0);
663 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
664 timer.currentStepLength(),
669 report += model_.wellModel().getDomainWellConvergence(domain, B_avg);
674 std::vector<int> getSubdomainOrder()
676 const auto& ebosSimulator = model_.ebosSimulator();
677 const auto& solution = ebosSimulator.model().solution(0);
679 std::vector<int> domain_order(domains_.size());
680 switch (model_.param().local_solve_approach_) {
681 case DomainSolveApproach::GaussSeidel: {
682 switch (model_.param().local_domain_ordering_) {
683 case DomainOrderingMeasure::AveragePressure: {
685 std::vector<std::pair<double, int>> avgpress_per_domain(domains_.size());
686 for (
const auto& domain : domains_) {
687 double press_sum = 0.0;
688 for (
const int c : domain.cells) {
689 press_sum += solution[c][Indices::pressureSwitchIdx];
691 const double avgpress = press_sum / domain.cells.size();
692 avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index);
695 std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end());
697 std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end());
698 for (std::size_t ii = 0; ii < domains_.size(); ++ii) {
699 domain_order[ii] = avgpress_per_domain[ii].second;
703 case DomainOrderingMeasure::Residual: {
705 const auto& residual = ebosSimulator.model().linearizer().residual();
706 const int num_vars = residual[0].size();
707 std::vector<std::pair<double, int>> maxres_per_domain(domains_.size());
708 for (
const auto& domain : domains_) {
710 for (
const int c : domain.cells) {
711 for (
int ii = 0; ii < num_vars; ++ii) {
712 maxres = std::max(maxres, std::fabs(residual[c][ii]));
715 maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index);
718 std::sort(maxres_per_domain.begin(), maxres_per_domain.end());
720 std::reverse(maxres_per_domain.begin(), maxres_per_domain.end());
721 for (std::size_t ii = 0; ii < domains_.size(); ++ii) {
722 domain_order[ii] = maxres_per_domain[ii].second;
730 case DomainSolveApproach::Jacobi:
732 std::iota(domain_order.begin(), domain_order.end(), 0);
739 template<
class GlobalEqVector>
740 void solveDomainJacobi(GlobalEqVector& solution,
741 GlobalEqVector& locally_solved,
742 SimulatorReportSingle& local_report,
744 const SimulatorTimerInterface& timer,
745 const Domain& domain)
747 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
748 auto initial_local_solution = Details::extractVector(solution, domain.cells);
749 auto res = solveDomain(domain, timer, iteration);
750 local_report = res.first;
751 if (local_report.converged) {
752 auto local_solution = Details::extractVector(solution, domain.cells);
753 Details::setGlobal(local_solution, domain.cells, locally_solved);
754 Details::setGlobal(initial_local_solution, domain.cells, solution);
755 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
757 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
758 Details::setGlobal(initial_local_solution, domain.cells, solution);
759 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
763 template<
class GlobalEqVector>
764 void solveDomainGaussSeidel(GlobalEqVector& solution,
765 GlobalEqVector& locally_solved,
766 SimulatorReportSingle& local_report,
768 const SimulatorTimerInterface& timer,
769 const Domain& domain)
771 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
772 auto initial_local_solution = Details::extractVector(solution, domain.cells);
773 auto res = solveDomain(domain, timer, iteration);
774 local_report = res.first;
775 if (!local_report.converged) {
778 const auto& convrep = res.second;
780 if (!convrep.wellFailed()) {
783 double cnv_sum = 0.0;
784 for (
const auto& rc : convrep.reservoirConvergence()) {
785 if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
786 mb_sum += rc.value();
787 }
else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) {
788 cnv_sum += rc.value();
792 const double acceptable_local_mb_sum = 1e-3;
793 const double acceptable_local_cnv_sum = 1.0;
794 if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) {
795 local_report.converged =
true;
796 OpmLog::debug(
"Accepting solution in unconverged domain " + std::to_string(domain.index));
800 if (local_report.converged) {
801 auto local_solution = Details::extractVector(solution, domain.cells);
802 Details::setGlobal(local_solution, domain.cells, locally_solved);
804 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
805 Details::setGlobal(initial_local_solution, domain.cells, solution);
806 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
810 double computeCnvErrorPvLocal(
const Domain& domain,
811 const std::vector<Scalar>& B_avg,
double dt)
const
814 const auto& ebosSimulator = model_.ebosSimulator();
815 const auto& ebosModel = ebosSimulator.model();
816 const auto& ebosProblem = ebosSimulator.problem();
817 const auto& ebosResid = ebosSimulator.model().linearizer().residual();
819 for (
const int cell_idx : domain.cells) {
820 const double pvValue = ebosProblem.referencePorosity(cell_idx, 0) *
821 ebosModel.dofTotalVolume(cell_idx);
822 const auto& cellResidual = ebosResid[cell_idx];
823 bool cnvViolated =
false;
825 for (
unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) {
827 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
828 cnvViolated = cnvViolated || (fabs(CNV) > model_.param().tolerance_cnv_);
838 BlackoilModelEbos<TypeTag>& model_;
839 std::vector<Domain> domains_;
840 std::vector<std::unique_ptr<Mat>> domain_matrices_;
841 std::vector<ISTLSolverType> domain_linsolvers_;
842 SimulatorReportSingle local_reports_accumulated_;