89 using BVector =
typename BlackoilModel<TypeTag>::BVector;
92 using Mat =
typename BlackoilModel<TypeTag>::Mat;
94 static constexpr int numEq = Indices::numEq;
102 , wellModel_(model.wellModel())
103 , rank_(model_.simulator().vanguard().grid().comm().rank())
127 model.
wellModel().setNlddAdapter(&wellModel_);
130 std::vector<int>
sizes(num_domains, 0);
136 using EntitySeed =
typename Grid::template Codim<0>::EntitySeed;
137 std::vector<std::vector<EntitySeed>>
seeds(num_domains);
138 std::vector<std::vector<int>>
partitions(num_domains);
146 const auto& grid = model_.simulator().vanguard().grid();
148 std::vector<int> count(num_domains, 0);
149 const auto& gridView = grid.leafGridView();
153 for (
auto it =
beg; it != end; ++it, ++cell) {
155 seeds[
p][count[
p]] = it->seed();
162 for (
int index = 0; index < num_domains; ++index) {
168 Dune::SubGridPart<Grid> view{grid, std::move(
seeds[index])};
172 this->domains_.emplace_back(index,
180 domain_matrices_.resize(num_domains);
183 for (
int index = 0; index < num_domains; ++index) {
187 const auto& eclState = model_.simulator().vanguard().eclState();
190 loc_param.init(eclState.getSimulationConfig().useCPR());
192 if (domains_[index].cells.size() < 200) {
195 loc_param.linear_solver_print_json_definition_ =
false;
198 domain_linsolvers_.back().setDomainIndex(index);
201 assert(
int(domains_.size()) == num_domains);
203 domain_reports_accumulated_.resize(num_domains);
209 local_reports_accumulated_,
210 domain_reports_accumulated_,
212 wellModel_.numLocalWellsEnd());
219 wellModel_.setupDomains(domains_);
223 template <
class NonlinearSolverType>
233 if (report.converged) {
242 auto& solution = model_.simulator().model().solution(0);
248 local_reports_accumulated_.success.pre_post_time +=
detailTimer.stop();
252 std::vector<SimulatorReportSingle>
domain_reports(domains_.size());
264 switch (model_.param().local_solve_approach_) {
265 case DomainSolveApproach::Jacobi:
270 case DomainSolveApproach::GaussSeidel:
285 logger.debug(fmt::format(
"Convergence failure in domain {} on rank {}." ,
domain.index, rank_));
303 int& num_domains =
counts[3];
311 if (
dr.total_newton_iterations == 0) {
317 domain_reports_accumulated_[i] +=
dr;
319 local_reports_accumulated_ +=
dr;
324 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
326 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0);
337 const auto& comm = model_.simulator().vanguard().grid().comm();
338 if (comm.size() > 1) {
339 const auto*
ccomm = model_.simulator().model().newtonMethod().linearSolver().comm();
342 ccomm->copyOwnerToAll(solution, solution);
345 const std::size_t
num = solution.size();
347 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
351 for (std::size_t
ii = 0;
ii <
num; ++
ii) {
356 model_.simulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(0);
365 OpmLog::debug(fmt::format(
"Local solves finished. Converged for {}/{} domains. {} domains did no work. {} total local Newton iterations.\n",
371 local_reports_accumulated_.success.pre_post_time +=
detailTimer.stop();
381 report.converged =
true;
389 return local_reports_accumulated_;
397 const auto dr_size = domain_reports_accumulated_.size();
400 domain_reports_accumulated_[i].success.num_wells = 0;
403 for (
const auto& [
wname,
domain] : wellModel_.well_domain()) {
404 domain_reports_accumulated_[
domain].success.num_wells++;
406 return domain_reports_accumulated_;
413 const auto& grid = this->model_.simulator().vanguard().grid();
414 const auto& elementMapper = this->model_.simulator().model().elementMapper();
415 const auto&
cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
428 const auto& grid = this->model_.simulator().vanguard().grid();
429 const auto& elementMapper = this->model_.simulator().model().elementMapper();
430 const auto&
cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
435 domain_reports_accumulated_,
444 solveDomain(
const Domain&
domain,
466 wellModel_.assemble(
modelSimulator.model().newtonMethod().numIterations(),
475 this->assembleReservoirDomain(
domain);
494 model_.wellModel().linearizeDomain(
domain,
502 const int max_iter = model_.param().max_local_solve_iterations_;
511 const int nc = grid.size(0);
515 this->solveJacobianSystemDomain(
domain, x);
516 model_.wellModel().postSolveDomain(x,
domain);
521 local_report.linear_solve_setup_time += model_.linearSolveSetupTime();
522 local_report.total_linear_iterations = model_.linearIterationsLastSolve();
527 this->updateDomainSolution(
domain, x);
538 wellModel_.assemble(
modelSimulator.model().newtonMethod().numIterations(),
546 this->assembleReservoirDomain(
domain);
561 model_.wellModel().linearizeDomain(
domain,
592 void assembleReservoirDomain(
const Domain&
domain)
596 model_.simulator().model().linearizer().linearizeDomain(
domain);
600 void solveJacobianSystemDomain(
const Domain&
domain, BVector&
global_x)
608 if (domain_matrices_[
domain.index]) {
613 auto&
jac = *domain_matrices_[
domain.index];
614 auto res = Details::extractVector(
modelSimulator.model().linearizer().residual(),
625 model_.linearSolveSetupTime() =
perfTimer.stop();
633 void updateDomainSolution(
const Domain&
domain,
const BVector&
dx)
636 auto& simulator = model_.simulator();
637 auto& newtonMethod = simulator.model().newtonMethod();
638 SolutionVector& solution = simulator.model().solution(0);
640 newtonMethod.update_(solution,
649 simulator.model().invalidateAndUpdateIntensiveQuantities(0,
domain);
653 std::pair<Scalar, Scalar> localDomainConvergenceData(
const Domain&
domain,
654 std::vector<Scalar>&
R_sum,
656 std::vector<Scalar>&
B_avg,
669 const auto& gridView =
domain.view;
670 const auto&
elemEndIt = gridView.template end<0>();
677 if (
elemIt->partitionType() != Dune::InteriorEntity) {
681 elemCtx.updatePrimaryStencil(
elem);
682 elemCtx.updatePrimaryIntensiveQuantities(0);
684 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
685 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
686 const auto& fs = intQuants.fluidState();
703 for (
int i = 0; i<
bSize; ++i )
711 ConvergenceReport getDomainReservoirConvergence(
const double reportTime,
716 std::vector<Scalar>&
B_avg,
719 using Vector = std::vector<Scalar>;
735 iteration >= model_.param().min_strict_cnv_iter_;
738 const Scalar
tol_cnv = model_.param().local_tolerance_scaling_cnv_
740 : model_.param().tolerance_cnv_);
743 const Scalar
tol_mb = model_.param().local_tolerance_scaling_mb_
744 * (
use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
757 ConvergenceReport report{reportTime};
758 using CR = ConvergenceReport;
761 CR::ReservoirFailure::Type
types[2] = { CR::ReservoirFailure::Type::MassBalance,
762 CR::ReservoirFailure::Type::Cnv };
764 for (
int ii : {0, 1}) {
765 if (std::isnan(
res[
ii])) {
766 report.setReservoirFailed({
types[
ii], CR::Severity::NotANumber,
compIdx});
767 logger.debug(
"NaN residual for " + model_.compNames().name(
compIdx) +
" equation.");
768 }
else if (
res[
ii] > model_.param().max_residual_allowed_) {
769 report.setReservoirFailed({
types[
ii], CR::Severity::TooLarge,
compIdx});
770 logger.debug(
"Too large residual for " + model_.compNames().name(
compIdx) +
" equation.");
771 }
else if (
res[
ii] < 0.0) {
772 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
773 logger.debug(
"Negative residual for " + model_.compNames().name(
compIdx) +
" equation.");
775 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
787 std::string
msg = fmt::format(
"Domain {} on rank {}, size {}, containing cell {}\n| Iter",
802 std::ostringstream
ss;
804 const std::streamsize
oprec =
ss.precision(3);
805 const std::ios::fmtflags
oflags =
ss.setf(std::ios::scientific);
821 ConvergenceReport getDomainConvergence(
const Domain&
domain,
822 const SimulatorTimerInterface& timer,
828 std::vector<Scalar>
B_avg(numEq, 0.0);
829 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
830 timer.currentStepLength(),
841 std::vector<int> getSubdomainOrder()
849 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
852 }
else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) {
855 switch (model_.param().local_domains_ordering_) {
856 case DomainOrderingMeasure::AveragePressure: {
858 for (
const auto&
domain : domains_) {
860 for (
const int c :
domain.cells) {
861 press_sum += solution[
c][Indices::pressureSwitchIdx];
868 case DomainOrderingMeasure::MaxPressure: {
870 for (
const auto&
domain : domains_) {
872 for (
const int c :
domain.cells) {
879 case DomainOrderingMeasure::Residual: {
881 const auto& residual =
modelSimulator.model().linearizer().residual();
882 const int num_vars = residual[0].size();
883 for (
const auto&
domain : domains_) {
885 for (
const int c :
domain.cells) {
899 [&
m](
const int i1,
const int i2){ return m[i1] > m[i2]; });
902 throw std::logic_error(
"Domain solve approach must be Jacobi or Gauss-Seidel");
906 template<
class GlobalEqVector>
907 void solveDomainJacobi(GlobalEqVector& solution,
912 const SimulatorTimerInterface& timer,
922 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
926 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
930 template<
class GlobalEqVector>
931 void solveDomainGaussSeidel(GlobalEqVector& solution,
936 const SimulatorTimerInterface& timer,
950 for (
const auto&
rc :
convrep.reservoirConvergence()) {
951 if (
rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
953 }
else if (
rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) {
963 logger.debug(fmt::format(
"Accepting solution in unconverged domain {} on rank {}.",
domain.index, rank_));
966 logger.debug(
"Unconverged local solution.");
969 logger.debug(
"Unconverged local solution with well convergence failures:");
970 for (
const auto&
wf :
convrep.wellFailures()) {
983 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0,
domain);
987 Scalar computeCnvErrorPvLocal(
const Domain&
domain,
988 const std::vector<Scalar>&
B_avg,
double dt)
const
991 const auto& simulator = model_.simulator();
992 const auto& model = simulator.model();
993 const auto& problem = simulator.problem();
994 const auto& residual = simulator.model().linearizer().residual();
1015 decltype(
auto) partitionCells()
const
1017 const auto& grid = this->model_.simulator().vanguard().grid();
1019 using GridView = std::remove_cv_t<std::remove_reference_t<
1020 decltype(grid.leafGridView())>>;
1022 using Element = std::remove_cv_t<std::remove_reference_t<
1023 typename GridView::template Codim<0>::Entity>>;
1025 const auto& param = this->model_.param();
1029 zoltan_ctrl.domain_imbalance = param.local_domains_partition_imbalance_;
1032 [elementMapper = &this->model_.simulator().model().elementMapper()]
1033 (
const Element& element)
1035 return elementMapper->index(element);
1039 [
cartMapper = &this->model_.simulator().vanguard().cartesianIndexMapper()]
1046 const auto need_wells = param.local_domains_partition_method_ ==
"zoltan";
1049 ? this->model_.simulator().vanguard().schedule().getWellsatEnd()
1050 : std::vector<Well>{};
1053 ? this->model_.simulator().vanguard().schedule().getPossibleFutureConnections()
1054 : std::unordered_map<std::string, std::set<int>> {};
1058 const int num_domains = (param.num_local_domains_ > 0)
1059 ? param.num_local_domains_
1061 return ::Opm::partitionCells(param.local_domains_partition_method_,
1062 num_domains, grid.leafGridView(), wells,
1064 param.local_domains_partition_well_neighbor_levels_);
1069 std::vector<Domain> domains_;
1070 std::vector<std::unique_ptr<Mat>> domain_matrices_;
1071 std::vector<ISTLSolverType> domain_linsolvers_;
1072 SimulatorReport local_reports_accumulated_;
1074 mutable std::vector<SimulatorReport> domain_reports_accumulated_;
Representing a part of a grid, in a way suitable for performing local solves.
Definition SubDomain.hpp:85