93 using BVector =
typename BlackoilModel<TypeTag>::BVector;
96 using Mat =
typename BlackoilModel<TypeTag>::Mat;
98 static constexpr int numEq = Indices::numEq;
106 , wellModel_(model.wellModel())
107 , rank_(model_.simulator().vanguard().grid().comm().rank())
131 model.
wellModel().setNlddAdapter(&wellModel_);
140 using EntitySeed =
typename Grid::template Codim<0>::EntitySeed;
150 const auto& grid = model_.simulator().vanguard().grid();
153 const auto& gridView = grid.leafGridView();
157 for (
auto it =
beg; it != end; ++it, ++cell) {
159 seeds[
p][count[
p]] = it->seed();
166 for (
int index = 0; index <
num_domains; ++index) {
172 Dune::SubGridPart<Grid> view{grid, std::move(
seeds[index])};
176 this->domains_.emplace_back(index,
187 for (
int index = 0; index <
num_domains; ++index) {
191 const auto& eclState = model_.simulator().vanguard().eclState();
194 loc_param.init(eclState.getSimulationConfig().useCPR());
196 if (domains_[index].cells.size() < 200) {
199 loc_param.linear_solver_print_json_definition_ =
false;
202 domain_linsolvers_.back().setDomainIndex(index);
212 wellModel_.setupDomains(domains_);
216 template <
class NonlinearSolverType>
227 if (report.converged) {
233 auto& solution = model_.simulator().model().solution(0);
242 std::vector<SimulatorReportSingle>
domain_reports(domains_.size());
252 switch (model_.param().local_solve_approach_) {
253 case DomainSolveApproach::Jacobi:
258 case DomainSolveApproach::GaussSeidel:
273 logger.debug(fmt::format(
"Convergence failure in domain {} on rank {}." ,
domain.index, rank_));
295 if (
dr.total_newton_iterations == 0) {
302 local_reports_accumulated_ +=
rep;
305 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
307 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0);
318 const auto& comm = model_.simulator().vanguard().grid().comm();
319 if (comm.size() > 1) {
320 const auto*
ccomm = model_.simulator().model().newtonMethod().linearSolver().comm();
323 ccomm->copyOwnerToAll(solution, solution);
326 const std::size_t
num = solution.size();
328 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
332 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
337 model_.simulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(0);
346 OpmLog::debug(fmt::format(
"Local solves finished. Converged for {}/{} domains. {} domains did no work. {} total local Newton iterations.\n",
358 report.converged =
true;
366 return local_reports_accumulated_;
373 const auto& elementMapper = this->model_.simulator().model().elementMapper();
374 const auto&
cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
378 const auto& grid = this->model_.simulator().vanguard().grid();
379 const auto& comm = grid.comm();
385 const auto p = this->reconstitutePartitionVector();
387 for (
const auto& cell :
elements(grid.leafGridView(), Dune::Partitions::interior)) {
395 if (this->rank_ == 0) {
396 auto fname =
odir /
"ResInsight_compatible_partition.txt";
409 const auto nDigit = 1 +
static_cast<int>(std::floor(std::log10(comm.size())));
414 for (
const auto& cell :
elements(grid.leafGridView(), Dune::Partitions::interior)) {
415 pfile << comm.rank() <<
' '
416 <<
cartMapper.cartesianIndex(elementMapper.index(cell)) <<
' '
417 <<
p[cell_index++] <<
'\n';
423 std::pair<SimulatorReportSingle, ConvergenceReport>
424 solveDomain(
const Domain&
domain,
433 Dune::Timer solveTimer;
449 wellModel_.assemble(
modelSimulator.model().newtonMethod().numIterations(),
453 report += this->assembleReservoirDomain(
domain);
462 report.converged =
true;
470 model_.wellModel().linearizeDomain(
domain,
474 report.assemble_time +=
tt1;
475 report.assemble_time_well +=
tt1;
478 const int max_iter = model_.param().max_local_solve_iterations_;
487 const int nc = grid.size(0);
491 this->solveJacobianSystemDomain(
domain, x);
492 model_.wellModel().postSolveDomain(x,
domain);
497 report.linear_solve_setup_time += model_.linearSolveSetupTime();
498 report.total_linear_iterations = model_.linearIterationsLastSolve();
503 this->updateDomainSolution(
domain, x);
514 wellModel_.assemble(
modelSimulator.model().newtonMethod().numIterations(),
517 report += this->assembleReservoirDomain(
domain);
531 model_.wellModel().linearizeDomain(
domain,
535 report.assemble_time +=
tt2;
536 report.assemble_time_well +=
tt2;
555 report.total_newton_iterations =
iter;
556 report.total_linearizations =
iter;
557 report.total_time = solveTimer.stop();
563 SimulatorReportSingle assembleReservoirDomain(
const Domain&
domain)
566 model_.simulator().model().linearizer().linearizeDomain(
domain);
567 return model_.wellModel().lastReport();
571 void solveJacobianSystemDomain(
const Domain&
domain, BVector&
global_x)
579 if (domain_matrices_[
domain.index]) {
584 auto&
jac = *domain_matrices_[
domain.index];
585 auto res = Details::extractVector(
modelSimulator.model().linearizer().residual(),
596 model_.linearSolveSetupTime() =
perfTimer.stop();
604 void updateDomainSolution(
const Domain&
domain,
const BVector&
dx)
606 auto& simulator = model_.simulator();
607 auto& newtonMethod = simulator.model().newtonMethod();
608 SolutionVector& solution = simulator.model().solution(0);
610 newtonMethod.update_(solution,
619 simulator.model().invalidateAndUpdateIntensiveQuantities(0,
domain);
623 std::pair<Scalar, Scalar> localDomainConvergenceData(
const Domain&
domain,
624 std::vector<Scalar>&
R_sum,
626 std::vector<Scalar>&
B_avg,
639 const auto& gridView =
domain.view;
640 const auto&
elemEndIt = gridView.template end<0>();
647 if (
elemIt->partitionType() != Dune::InteriorEntity) {
651 elemCtx.updatePrimaryStencil(
elem);
652 elemCtx.updatePrimaryIntensiveQuantities(0);
654 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
655 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
656 const auto& fs = intQuants.fluidState();
673 for (
int i = 0; i<
bSize; ++i )
681 ConvergenceReport getDomainReservoirConvergence(
const double reportTime,
686 std::vector<Scalar>&
B_avg,
689 using Vector = std::vector<Scalar>;
705 iteration >= model_.param().min_strict_cnv_iter_;
708 const Scalar
tol_cnv = model_.param().local_tolerance_scaling_cnv_
710 : model_.param().tolerance_cnv_);
713 const Scalar
tol_mb = model_.param().local_tolerance_scaling_mb_
714 * (
use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
727 ConvergenceReport report{reportTime};
728 using CR = ConvergenceReport;
731 CR::ReservoirFailure::Type
types[2] = { CR::ReservoirFailure::Type::MassBalance,
732 CR::ReservoirFailure::Type::Cnv };
734 for (
int ii : {0, 1}) {
735 if (std::isnan(
res[
ii])) {
736 report.setReservoirFailed({
types[
ii], CR::Severity::NotANumber,
compIdx});
737 logger.debug(
"NaN residual for " + model_.compNames().name(
compIdx) +
" equation.");
738 }
else if (
res[
ii] > model_.param().max_residual_allowed_) {
739 report.setReservoirFailed({
types[
ii], CR::Severity::TooLarge,
compIdx});
740 logger.debug(
"Too large residual for " + model_.compNames().name(
compIdx) +
" equation.");
741 }
else if (
res[
ii] < 0.0) {
742 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
743 logger.debug(
"Negative residual for " + model_.compNames().name(
compIdx) +
" equation.");
745 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
757 std::string
msg = fmt::format(
"Domain {} on rank {}, size {}, containing cell {}\n| Iter",
772 std::ostringstream
ss;
774 const std::streamsize
oprec =
ss.precision(3);
775 const std::ios::fmtflags
oflags =
ss.setf(std::ios::scientific);
791 ConvergenceReport getDomainConvergence(
const Domain&
domain,
792 const SimulatorTimerInterface& timer,
797 std::vector<Scalar>
B_avg(numEq, 0.0);
798 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
799 timer.currentStepLength(),
810 std::vector<int> getSubdomainOrder()
818 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
821 }
else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) {
824 switch (model_.param().local_domains_ordering_) {
825 case DomainOrderingMeasure::AveragePressure: {
827 for (
const auto&
domain : domains_) {
829 for (
const int c :
domain.cells) {
830 press_sum += solution[
c][Indices::pressureSwitchIdx];
837 case DomainOrderingMeasure::MaxPressure: {
839 for (
const auto&
domain : domains_) {
841 for (
const int c :
domain.cells) {
848 case DomainOrderingMeasure::Residual: {
850 const auto& residual =
modelSimulator.model().linearizer().residual();
851 const int num_vars = residual[0].size();
852 for (
const auto&
domain : domains_) {
854 for (
const int c :
domain.cells) {
868 [&
m](
const int i1,
const int i2){ return m[i1] > m[i2]; });
871 throw std::logic_error(
"Domain solve approach must be Jacobi or Gauss-Seidel");
875 template<
class GlobalEqVector>
876 void solveDomainJacobi(GlobalEqVector& solution,
881 const SimulatorTimerInterface& timer,
892 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
896 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
900 template<
class GlobalEqVector>
901 void solveDomainGaussSeidel(GlobalEqVector& solution,
906 const SimulatorTimerInterface& timer,
922 for (
const auto&
rc :
convrep.reservoirConvergence()) {
923 if (
rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
925 }
else if (
rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) {
934 logger.debug(fmt::format(
"Accepting solution in unconverged domain {} on rank {}.",
domain.index, rank_));
937 logger.debug(
"Unconverged local solution.");
940 logger.debug(
"Unconverged local solution with well convergence failures:");
941 for (
const auto&
wf :
convrep.wellFailures()) {
952 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
956 Scalar computeCnvErrorPvLocal(
const Domain&
domain,
957 const std::vector<Scalar>&
B_avg,
double dt)
const
960 const auto& simulator = model_.simulator();
961 const auto& model = simulator.model();
962 const auto& problem = simulator.problem();
963 const auto& residual = simulator.model().linearizer().residual();
984 decltype(
auto) partitionCells()
const
986 const auto& grid = this->model_.simulator().vanguard().grid();
988 using GridView = std::remove_cv_t<std::remove_reference_t<
989 decltype(grid.leafGridView())>>;
991 using Element = std::remove_cv_t<std::remove_reference_t<
992 typename GridView::template Codim<0>::Entity>>;
994 const auto& param = this->model_.param();
998 zoltan_ctrl.domain_imbalance = param.local_domains_partition_imbalance_;
1001 [elementMapper = &this->model_.simulator().model().elementMapper()]
1002 (
const Element& element)
1004 return elementMapper->index(element);
1008 [
cartMapper = &this->model_.simulator().vanguard().cartesianIndexMapper()]
1015 const auto need_wells = param.local_domains_partition_method_ ==
"zoltan";
1018 ? this->model_.simulator().vanguard().schedule().getWellsatEnd()
1019 : std::vector<Well>{};
1022 ? this->model_.simulator().vanguard().schedule().getPossibleFutureConnections()
1023 : std::unordered_map<std::string, std::set<int>> {};
1027 const int num_domains = (param.num_local_domains_ > 0)
1028 ? param.num_local_domains_
1030 return ::Opm::partitionCells(param.local_domains_partition_method_,
1033 param.local_domains_partition_well_neighbor_levels_);
1036 std::vector<int> reconstitutePartitionVector()
const
1038 const auto& grid = this->model_.simulator().vanguard().grid();
1040 auto numD = std::vector<int>(grid.comm().size() + 1, 0);
1041 numD[grid.comm().rank() + 1] =
static_cast<int>(this->domains_.size());
1042 grid.comm().sum(
numD.data(),
numD.size());
1043 std::partial_sum(
numD.begin(),
numD.end(),
numD.begin());
1045 auto p = std::vector<int>(grid.size(0));
1046 auto maxCellIdx = std::numeric_limits<int>::min();
1048 auto d =
numD[grid.comm().rank()];
1050 for (
const auto& cell :
domain.cells) {
1066 std::vector<Domain> domains_;
1067 std::vector<std::unique_ptr<Mat>> domain_matrices_;
1068 std::vector<ISTLSolverType> domain_linsolvers_;
1069 SimulatorReportSingle local_reports_accumulated_;