24#ifndef OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED
27#include <dune/common/timer.hh>
29#include <opm/grid/common/SubGridPart.hpp>
57#include <fmt/format.h>
76template<
class TypeTag>
class BlackoilModel;
79template <
class TypeTag>
96 static constexpr int numEq = Indices::numEq;
104 , wellModel_(model.wellModel())
105 , rank_(model_.simulator().vanguard().grid().comm().rank())
108 const auto& [partition_vector_initial, num_domains_initial] = this->partitionCells();
110 int num_domains = num_domains_initial;
111 std::vector<int> partition_vector = partition_vector_initial;
117 bool isolated_cells =
false;
118 for (
auto& domainId : partition_vector) {
120 domainId = num_domains;
121 isolated_cells =
true;
124 if (isolated_cells) {
129 model.
wellModel().setNlddAdapter(&wellModel_);
132 std::vector<int> sizes(num_domains, 0);
133 for (
const auto& p : partition_vector) {
138 using EntitySeed =
typename Grid::template Codim<0>::EntitySeed;
139 std::vector<std::vector<EntitySeed>> seeds(num_domains);
140 std::vector<std::vector<int>> partitions(num_domains);
141 for (
int domain = 0; domain < num_domains; ++domain) {
142 seeds[domain].resize(sizes[domain]);
143 partitions[domain].resize(sizes[domain]);
148 const auto& grid = model_.simulator().vanguard().grid();
150 std::vector<int> count(num_domains, 0);
151 const auto& gridView = grid.leafGridView();
152 const auto beg = gridView.template begin<0, Dune::Interior_Partition>();
153 const auto end = gridView.template end<0, Dune::Interior_Partition>();
155 for (
auto it = beg; it != end; ++it, ++cell) {
156 const int p = partition_vector[cell];
157 seeds[p][count[p]] = it->seed();
158 partitions[p][count[p]] = cell;
161 assert(count == sizes);
164 for (
int index = 0; index < num_domains; ++index) {
165 std::vector<bool> interior(partition_vector.size(),
false);
166 for (
int ix : partitions[index]) {
170 Dune::SubGridPart<Grid> view{grid, std::move(seeds[index])};
173 const bool skip = isolated_cells && (index == num_domains - 1);
174 this->domains_.emplace_back(index,
175 std::move(partitions[index]),
182 const auto numCells = grid.size(0);
183 previousMobilities_.resize(numCells * FluidSystem::numActivePhases(), 0.0);
184 for (
const auto& domain : domains_) {
185 updateMobilities(domain);
189 domain_needs_solving_.resize(num_domains,
true);
192 domain_matrices_.resize(num_domains);
195 for (
int index = 0; index < num_domains; ++index) {
199 const auto& eclState = model_.simulator().vanguard().eclState();
202 loc_param.
init(eclState.getSimulationConfig().useCPR());
204 if (domains_[index].cells.size() < 200) {
208 const bool force_serial =
true;
209 domain_linsolvers_.emplace_back(model_.simulator(), loc_param, force_serial);
210 domain_linsolvers_.back().setDomainIndex(index);
213 assert(
int(domains_.size()) == num_domains);
215 domain_reports_accumulated_.resize(num_domains);
221 local_reports_accumulated_,
222 domain_reports_accumulated_,
224 wellModel_.numLocalWellsEnd());
231 wellModel_.setupDomains(domains_);
235 template <
class NonlinearSolverType>
238 NonlinearSolverType& nonlinear_solver)
243 if (iteration < model_.param().nldd_num_initial_newton_iter_) {
244 report = model_.nonlinearIterationNewton(iteration,
250 model_.initialLinearization(report, iteration, model_.param().newton_min_iter_,
251 model_.param().newton_max_iter_, timer);
258 Dune::Timer localSolveTimer;
259 Dune::Timer detailTimer;
260 localSolveTimer.start();
262 auto& solution = model_.simulator().model().solution(0);
263 auto initial_solution = solution;
264 auto locally_solved = initial_solution;
267 const auto domain_order = this->getSubdomainOrder();
272 std::vector<SimulatorReportSingle> domain_reports(domains_.size());
273 for (
const int domain_index : domain_order) {
274 const auto& domain = domains_[domain_index];
279 domain_needs_solving_[domain_index] = checkIfSubdomainNeedsSolving(domain, iteration);
281 updateMobilities(domain);
283 if (domain.skip || !domain_needs_solving_[domain_index]) {
286 domain_reports[domain.index] = local_report;
290 switch (model_.param().local_solve_approach_) {
292 solveDomainJacobi(solution, locally_solved, local_report, logger,
293 iteration, timer, domain);
297 solveDomainGaussSeidel(solution, locally_solved, local_report, logger,
298 iteration, timer, domain);
311 logger.
debug(fmt::format(
"Convergence failure in domain {} on rank {}." , domain.index, rank_));
314 domain_reports[domain.index] = local_report;
320 global_logger.logMessages();
325 std::array<int, 5> counts{ 0, 0, 0,
static_cast<int>(domain_reports.size()), 0 };
326 int& num_converged = counts[0];
327 int& num_converged_already = counts[1];
328 int& num_local_newtons = counts[2];
329 int& num_domains = counts[3];
330 int& num_skipped = counts[4];
332 auto step_newtons = 0;
333 const auto dr_size = domain_reports.size();
334 for (
auto i = 0*dr_size; i < dr_size; ++i) {
336 domain_needs_solving_[i] =
false;
337 const auto& dr = domain_reports[i];
340 if (dr.total_newton_iterations == 0) {
341 ++num_converged_already;
345 domain_needs_solving_[i] =
true;
348 if (dr.skipped_domains) {
351 step_newtons += dr.total_newton_iterations;
353 domain_reports_accumulated_[i] += dr;
355 local_reports_accumulated_ += dr;
357 num_local_newtons = step_newtons;
361 solution = locally_solved;
362 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0);
373 const auto& comm = model_.simulator().vanguard().grid().comm();
374 if (comm.size() > 1) {
375 const auto* ccomm = model_.simulator().model().newtonMethod().linearSolver().comm();
378 ccomm->copyOwnerToAll(solution, solution);
381 const std::size_t num = solution.size();
382 Dune::BlockVector<std::size_t> allmeanings(num);
383 for (std::size_t ii = 0; ii < num; ++ii) {
386 ccomm->copyOwnerToAll(allmeanings, allmeanings);
387 for (std::size_t ii = 0; ii < num; ++ii) {
392 model_.simulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(0);
395 comm.sum(counts.data(), counts.size());
399 const bool is_iorank = this->rank_ == 0;
401 OpmLog::debug(fmt::format(
"Local solves finished. Converged for {}/{} domains. {} domains were skipped. {} domains did no work. {} total local Newton iterations.\n",
402 num_converged, num_domains, num_skipped, num_converged_already, num_local_newtons));
404 auto total_local_solve_time = localSolveTimer.stop();
414 auto rep = model_.nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver);
425 return local_reports_accumulated_;
433 const auto dr_size = domain_reports_accumulated_.size();
435 for (
auto i = 0*dr_size; i < dr_size; ++i) {
436 domain_reports_accumulated_[i].success.num_wells = 0;
439 for (
const auto& [wname, domain] : wellModel_.well_domain()) {
440 domain_reports_accumulated_[domain].success.num_wells++;
442 return domain_reports_accumulated_;
449 const auto& grid = this->model_.simulator().vanguard().grid();
450 const auto& elementMapper = this->model_.simulator().model().elementMapper();
451 const auto& cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
464 const auto& grid = this->model_.simulator().vanguard().grid();
465 const auto& elementMapper = this->model_.simulator().model().elementMapper();
466 const auto& cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
471 domain_reports_accumulated_,
480 solveDomain(
const Domain& domain,
484 [[maybe_unused]]
const int global_iteration,
485 const bool initial_assembly_required)
487 auto& modelSimulator = model_.simulator();
488 Dune::Timer detailTimer;
490 modelSimulator.model().newtonMethod().setIterationIndex(0);
497 if (initial_assembly_required) {
499 modelSimulator.model().newtonMethod().setIterationIndex(iter);
502 wellModel_.assemble(modelSimulator.model().newtonMethod().numIterations(),
503 modelSimulator.timeStepSize(),
505 const double tt0 = detailTimer.stop();
511 this->assembleReservoirDomain(domain);
517 std::vector<Scalar> resnorms;
518 auto convreport = this->getDomainConvergence(domain, timer, 0, logger, resnorms);
520 if (convreport.converged()) {
530 model_.wellModel().linearizeDomain(domain,
531 modelSimulator.model().linearizer().jacobian(),
532 modelSimulator.model().linearizer().residual());
533 const double tt1 = detailTimer.stop();
538 const int max_iter = model_.param().max_local_solve_iterations_;
539 const auto& grid = modelSimulator.vanguard().grid();
540 double damping_factor = 1.0;
541 std::vector<std::vector<Scalar>> convergence_history;
542 convergence_history.reserve(20);
543 convergence_history.push_back(resnorms);
547 const int nc = grid.size(0);
551 this->solveJacobianSystemDomain(domain, x);
552 model_.wellModel().postSolveDomain(x, domain);
553 if (damping_factor != 1.0) {
563 this->updateDomainSolution(domain, x);
570 modelSimulator.model().newtonMethod().setIterationIndex(iter);
574 wellModel_.assemble(modelSimulator.model().newtonMethod().numIterations(),
575 modelSimulator.timeStepSize(),
577 const double tt3 = detailTimer.stop();
582 this->assembleReservoirDomain(domain);
589 convreport = this->getDomainConvergence(domain, timer, iter, logger, resnorms);
590 convergence_history.push_back(resnorms);
597 model_.wellModel().linearizeDomain(domain,
598 modelSimulator.model().linearizer().jacobian(),
599 modelSimulator.model().linearizer().residual());
600 const double tt2 = detailTimer.stop();
605 if (!convreport.converged() && !convreport.wellFailed()) {
606 bool oscillate =
false;
607 bool stagnate =
false;
608 const auto num_residuals = convergence_history.front().size();
610 Scalar{0.2}, 1, oscillate, stagnate);
612 damping_factor *= 0.85;
613 logger.
debug(fmt::format(
"| Damping factor is now {}", damping_factor));
616 }
while (!convreport.converged() && iter <= max_iter);
618 modelSimulator.problem().endIteration();
620 local_report.
converged = convreport.converged();
628 void assembleReservoirDomain(
const Domain& domain)
630 OPM_TIMEBLOCK(assembleReservoirDomain);
632 model_.simulator().model().linearizer().linearizeDomain(domain);
636 void solveJacobianSystemDomain(
const Domain& domain,
BVector& global_x)
638 const auto& modelSimulator = model_.simulator();
640 Dune::Timer perfTimer;
643 const Mat& main_matrix = modelSimulator.model().linearizer().jacobian().istlMatrix();
644 if (domain_matrices_[domain.index]) {
647 domain_matrices_[domain.index] = std::make_unique<Mat>(
Details::extractMatrix(main_matrix, domain.cells));
649 auto& jac = *domain_matrices_[domain.index];
658 auto& linsolver = domain_linsolvers_[domain.index];
660 linsolver.prepare(jac, res);
661 model_.linearSolveSetupTime() = perfTimer.stop();
662 linsolver.setResidual(res);
669 void updateDomainSolution(
const Domain& domain,
const BVector& dx)
671 OPM_TIMEBLOCK(updateDomainSolution);
672 auto& simulator = model_.simulator();
673 auto& newtonMethod = simulator.model().newtonMethod();
676 newtonMethod.update_(solution,
685 simulator.model().invalidateAndUpdateIntensiveQuantities(0, domain);
689 std::pair<Scalar, Scalar> localDomainConvergenceData(
const Domain& domain,
690 std::vector<Scalar>& R_sum,
691 std::vector<Scalar>& maxCoeff,
692 std::vector<Scalar>& B_avg,
693 std::vector<int>& maxCoeffCell)
695 const auto& modelSimulator = model_.simulator();
698 Scalar numAquiferPvSumLocal = 0.0;
699 const auto& model = modelSimulator.model();
700 const auto& problem = modelSimulator.problem();
702 const auto& modelResid = modelSimulator.model().linearizer().residual();
705 const auto& gridView = domain.view;
706 const auto& elemEndIt = gridView.template end<0>();
707 IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
709 for (
auto elemIt = gridView.template begin</*codim=*/0>();
713 if (elemIt->partitionType() != Dune::InteriorEntity) {
716 const auto& elem = *elemIt;
717 elemCtx.updatePrimaryStencil(elem);
718 elemCtx.updatePrimaryIntensiveQuantities(0);
720 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
721 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
722 const auto& fs = intQuants.fluidState();
724 const auto pvValue = problem.referencePorosity(cell_idx, 0) *
725 model.dofTotalVolume(cell_idx);
726 pvSumLocal += pvValue;
728 if (isNumericalAquiferCell(elem))
730 numAquiferPvSumLocal += pvValue;
733 model_.getMaxCoeff(cell_idx, intQuants, fs, modelResid, pvValue,
734 B_avg, R_sum, maxCoeff, maxCoeffCell);
738 const int bSize = B_avg.size();
739 for (
int i = 0; i<bSize; ++i )
741 B_avg[ i ] /=
Scalar(domain.cells.size());
744 return {pvSumLocal, numAquiferPvSumLocal};
747 ConvergenceReport getDomainReservoirConvergence(
const double reportTime,
751 DeferredLogger& logger,
752 std::vector<Scalar>& B_avg,
753 std::vector<Scalar>& residual_norms)
755 using Vector = std::vector<Scalar>;
757 const int numComp =
numEq;
758 Vector R_sum(numComp, 0.0 );
759 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest() );
760 std::vector<int> maxCoeffCell(numComp, -1);
761 const auto [ pvSum, numAquiferPvSum]
762 = this->localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell);
764 auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt);
765 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
770 const bool use_relaxed_cnv = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
771 iteration >= model_.param().min_strict_cnv_iter_;
774 const Scalar tol_cnv = model_.param().local_tolerance_scaling_cnv_
775 * (use_relaxed_cnv ? model_.param().tolerance_cnv_relaxed_
776 : model_.param().tolerance_cnv_);
778 const bool use_relaxed_mb = iteration >= model_.param().min_strict_mb_iter_;
779 const Scalar tol_mb = model_.param().local_tolerance_scaling_mb_
780 * (use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
783 std::vector<Scalar> CNV(numComp);
784 std::vector<Scalar> mass_balance_residual(numComp);
785 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
787 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
788 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
789 residual_norms.push_back(CNV[compIdx]);
793 ConvergenceReport report{reportTime};
794 using CR = ConvergenceReport;
795 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
796 Scalar res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
797 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
798 CR::ReservoirFailure::Type::Cnv };
799 Scalar tol[2] = { tol_mb, tol_cnv };
800 for (
int ii : {0, 1}) {
801 if (std::isnan(res[ii])) {
802 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
803 logger.debug(
"NaN residual for " + model_.compNames().name(compIdx) +
" equation.");
804 }
else if (res[ii] > model_.param().max_residual_allowed_) {
805 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
806 logger.debug(
"Too large residual for " + model_.compNames().name(compIdx) +
" equation.");
807 }
else if (res[ii] < 0.0) {
808 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
809 logger.debug(
"Negative residual for " + model_.compNames().name(compIdx) +
" equation.");
810 }
else if (res[ii] > tol[ii]) {
811 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
814 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
819 const bool converged_at_initial_state = (report.converged() && iteration == 0);
820 if (!converged_at_initial_state) {
821 if (iteration == 0) {
823 std::string msg = fmt::format(
"Domain {} on rank {}, size {}, containing cell {}\n| Iter",
824 domain.index, this->rank_, domain.cells.size(), domain.cells[0]);
825 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
827 msg += model_.compNames().name(compIdx)[0];
830 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
832 msg += model_.compNames().name(compIdx)[0];
838 std::ostringstream ss;
840 const std::streamsize oprec = ss.precision(3);
841 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
842 ss << std::setw(4) << iteration;
843 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
844 ss << std::setw(11) << mass_balance_residual[compIdx];
846 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
847 ss << std::setw(11) << CNV[compIdx];
851 logger.debug(ss.str());
857 ConvergenceReport getDomainConvergence(
const Domain& domain,
858 const SimulatorTimerInterface& timer,
860 DeferredLogger& logger,
861 std::vector<Scalar>& residual_norms)
863 OPM_TIMEBLOCK(getDomainConvergence);
864 std::vector<Scalar> B_avg(
numEq, 0.0);
865 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
866 timer.currentStepLength(),
872 report += wellModel_.getWellConvergence(domain, B_avg, logger);
877 std::vector<int> getSubdomainOrder()
879 const auto& modelSimulator = model_.simulator();
880 const auto& solution = modelSimulator.model().solution(0);
882 std::vector<int> domain_order(domains_.size());
883 std::iota(domain_order.begin(), domain_order.end(), 0);
890 std::vector<Scalar> measure_per_domain(domains_.size());
891 switch (model_.param().local_domains_ordering_) {
894 for (
const auto& domain : domains_) {
896 std::accumulate(domain.cells.begin(), domain.cells.end(),
Scalar{0},
897 [&solution](
const auto acc,
const auto c)
898 { return acc + solution[c][Indices::pressureSwitchIdx]; });
899 const Scalar avgpress = press_sum / domain.cells.size();
900 measure_per_domain[domain.index] = avgpress;
906 for (
const auto& domain : domains_) {
907 measure_per_domain[domain.index] =
908 std::accumulate(domain.cells.begin(), domain.cells.end(),
Scalar{0},
909 [&solution](
const auto acc,
const auto c)
910 { return std::max(acc, solution[c][Indices::pressureSwitchIdx]); });
916 const auto& residual = modelSimulator.model().linearizer().residual();
917 const int num_vars = residual[0].size();
918 for (
const auto& domain : domains_) {
920 for (
const int c : domain.cells) {
921 for (
int ii = 0; ii < num_vars; ++ii) {
922 maxres = std::max(maxres, std::fabs(residual[c][ii]));
925 measure_per_domain[domain.index] = maxres;
932 const auto& m = measure_per_domain;
933 std::stable_sort(domain_order.begin(), domain_order.end(),
934 [&m](
const int i1,
const int i2){ return m[i1] > m[i2]; });
937 throw std::logic_error(
"Domain solve approach must be Jacobi or Gauss-Seidel");
941 template<
class GlobalEqVector>
942 void solveDomainJacobi(GlobalEqVector& solution,
943 GlobalEqVector& locally_solved,
944 SimulatorReportSingle& local_report,
945 DeferredLogger& logger,
947 const SimulatorTimerInterface& timer,
950 auto initial_local_well_primary_vars = wellModel_.getPrimaryVarsDomain(domain.index);
952 auto convrep = solveDomain(domain, timer, local_report, logger, iteration,
false);
953 if (local_report.converged) {
957 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
959 wellModel_.setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars);
961 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
965 template<
class GlobalEqVector>
966 void solveDomainGaussSeidel(GlobalEqVector& solution,
967 GlobalEqVector& locally_solved,
968 SimulatorReportSingle& local_report,
969 DeferredLogger& logger,
971 const SimulatorTimerInterface& timer,
974 auto initial_local_well_primary_vars = wellModel_.getPrimaryVarsDomain(domain.index);
976 auto convrep = solveDomain(domain, timer, local_report, logger, iteration,
true);
977 if (!local_report.converged) {
981 if (!convrep.wellFailed()) {
985 for (
const auto& rc : convrep.reservoirConvergence()) {
987 mb_sum += rc.value();
989 cnv_sum += rc.value();
993 const Scalar acceptable_local_mb_sum = 1e-3;
994 const Scalar acceptable_local_cnv_sum = 1.0;
995 if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) {
996 local_report.converged =
true;
997 local_report.accepted_unconverged_domains += 1;
998 logger.debug(fmt::format(
"Accepting solution in unconverged domain {} on rank {}.", domain.index, rank_));
999 logger.debug(fmt::format(
"Value of mb_sum: {} cnv_sum: {}", mb_sum, cnv_sum));
1001 logger.debug(
"Unconverged local solution.");
1004 logger.debug(
"Unconverged local solution with well convergence failures:");
1005 for (
const auto& wf : convrep.wellFailures()) {
1010 if (local_report.converged) {
1011 local_report.converged_domains += 1;
1015 local_report.unconverged_domains += 1;
1016 wellModel_.setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars);
1018 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(0, domain);
1023 const std::vector<Scalar>& B_avg,
double dt)
const
1026 const auto& simulator = model_.simulator();
1027 const auto& model = simulator.model();
1028 const auto& problem = simulator.problem();
1029 const auto& residual = simulator.model().linearizer().residual();
1031 for (
const int cell_idx : domain.cells) {
1032 const Scalar pvValue = problem.referencePorosity(cell_idx, 0) *
1033 model.dofTotalVolume(cell_idx);
1034 const auto& cellResidual = residual[cell_idx];
1035 bool cnvViolated =
false;
1037 for (
unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) {
1039 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
1040 cnvViolated = cnvViolated || (fabs(CNV) > model_.param().tolerance_cnv_);
1050 decltype(
auto) partitionCells()
const
1052 const auto& grid = this->model_.simulator().vanguard().grid();
1054 using GridView = std::remove_cv_t<std::remove_reference_t<
1055 decltype(grid.leafGridView())>>;
1057 using Element = std::remove_cv_t<std::remove_reference_t<
1058 typename GridView::template Codim<0>::Entity>>;
1060 const auto& param = this->model_.param();
1062 auto zoltan_ctrl = ZoltanPartitioningControl<Element>{};
1064 zoltan_ctrl.domain_imbalance = param.local_domains_partition_imbalance_;
1067 [elementMapper = &this->model_.simulator().model().elementMapper()]
1068 (
const Element& element)
1070 return elementMapper->index(element);
1073 zoltan_ctrl.local_to_global =
1074 [cartMapper = &this->model_.simulator().vanguard().cartesianIndexMapper()]
1077 return cartMapper->cartesianIndex(elemIdx);
1081 const auto need_wells = param.local_domains_partition_method_ ==
"zoltan";
1083 const auto wells = need_wells
1084 ? this->model_.simulator().vanguard().schedule().getWellsatEnd()
1085 : std::vector<Well>{};
1087 const auto& possibleFutureConnectionSet = need_wells
1088 ? this->model_.simulator().vanguard().schedule().getPossibleFutureConnections()
1089 : std::unordered_map<std::string, std::set<int>> {};
1092 constexpr int default_cells_per_domain = 1000;
1093 const int num_domains = (param.num_local_domains_ > 0)
1094 ? param.num_local_domains_
1097 num_domains, grid.leafGridView(), wells,
1098 possibleFutureConnectionSet, zoltan_ctrl,
1099 param.local_domains_partition_well_neighbor_levels_);
1102 void updateMobilities(
const Domain& domain)
1104 if (domain.skip || model_.param().nldd_relative_mobility_change_tol_ == 0.0) {
1107 const auto numActivePhases = FluidSystem::numActivePhases();
1108 for (
const auto globalDofIdx : domain.cells) {
1109 const auto& intQuants = model_.simulator().model().intensiveQuantities(globalDofIdx, 0);
1111 for (
unsigned activePhaseIdx = 0; activePhaseIdx < numActivePhases; ++activePhaseIdx) {
1112 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
1113 const auto mobIdx = globalDofIdx * numActivePhases + activePhaseIdx;
1114 previousMobilities_[mobIdx] = getValue(intQuants.mobility(phaseIdx));
1119 bool checkIfSubdomainNeedsSolving(
const Domain& domain,
const int iteration)
1127 if (domain_needs_solving_[domain.index]) {
1132 if (model_.param().nldd_relative_mobility_change_tol_ == 0.0) {
1137 if (iteration == 0) {
1141 return checkSubdomainChangeRelative(domain);
1144 bool checkSubdomainChangeRelative(
const Domain& domain)
1146 const auto numActivePhases = FluidSystem::numActivePhases();
1149 for (
const auto globalDofIdx : domain.cells) {
1150 const auto& intQuants = model_.simulator().model().intensiveQuantities(globalDofIdx, 0);
1154 for (
unsigned activePhaseIdx = 0; activePhaseIdx < numActivePhases; ++activePhaseIdx) {
1155 const auto mobIdx = globalDofIdx * numActivePhases + activePhaseIdx;
1156 cellMob += previousMobilities_[mobIdx] / numActivePhases;
1160 for (
unsigned activePhaseIdx = 0; activePhaseIdx < numActivePhases; ++activePhaseIdx) {
1161 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
1162 const auto mobIdx = globalDofIdx * numActivePhases + activePhaseIdx;
1163 const auto mobility = getValue(intQuants.mobility(phaseIdx));
1164 const auto relDiff = std::abs(mobility - previousMobilities_[mobIdx]) / cellMob;
1165 if (relDiff > model_.param().nldd_relative_mobility_change_tol_) {
1173 BlackoilModel<TypeTag>& model_;
1174 BlackoilWellModelNldd<TypeTag> wellModel_;
1175 std::vector<Domain> domains_;
1176 std::vector<std::unique_ptr<Mat>> domain_matrices_;
1177 std::vector<ISTLSolverType> domain_linsolvers_;
1178 SimulatorReport local_reports_accumulated_;
1180 mutable std::vector<SimulatorReport> domain_reports_accumulated_;
1183 std::vector<Scalar> previousMobilities_;
1185 std::vector<bool> domain_needs_solving_;
A NLDD implementation for three-phase black oil.
Definition: BlackoilModelNldd.hpp:81
GetPropType< TypeTag, Properties::SolutionVector > SolutionVector
Definition: BlackoilModelNldd.hpp:89
GetPropType< TypeTag, Properties::Indices > Indices
Definition: BlackoilModelNldd.hpp:86
const std::vector< SimulatorReport > & domainAccumulatedReports() const
return the statistics of local solves accumulated for each domain on this rank
Definition: BlackoilModelNldd.hpp:429
void writePartitions(const std::filesystem::path &odir) const
Definition: BlackoilModelNldd.hpp:447
BlackoilModelNldd(BlackoilModel< TypeTag > &model)
The constructor sets up the subdomains.
Definition: BlackoilModelNldd.hpp:102
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: BlackoilModelNldd.hpp:87
typename BlackoilModel< TypeTag >::Mat Mat
Definition: BlackoilModelNldd.hpp:94
const SimulatorReport & localAccumulatedReports() const
return the statistics of local solves accumulated for this rank
Definition: BlackoilModelNldd.hpp:423
SimulatorReportSingle nonlinearIterationNldd(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Do one non-linear NLDD iteration.
Definition: BlackoilModelNldd.hpp:236
SubDomain< Grid > Domain
Definition: BlackoilModelNldd.hpp:92
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: BlackoilModelNldd.hpp:83
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: BlackoilModelNldd.hpp:84
void prepareStep()
Called before starting a time step.
Definition: BlackoilModelNldd.hpp:228
void writeNonlinearIterationsPerCell(const std::filesystem::path &odir) const
Write the number of nonlinear iterations per cell to a file in ResInsight compatible format.
Definition: BlackoilModelNldd.hpp:462
static constexpr int numEq
Definition: BlackoilModelNldd.hpp:96
GetPropType< TypeTag, Properties::Grid > Grid
Definition: BlackoilModelNldd.hpp:85
typename BlackoilModel< TypeTag >::BVector BVector
Definition: BlackoilModelNldd.hpp:91
Definition: BlackoilModel.hpp:61
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModel.hpp:282
typename SparseMatrixAdapter::IstlMatrix Mat
Definition: BlackoilModel.hpp:106
Dune::BlockVector< VectorBlockType > BVector
Definition: BlackoilModel.hpp:107
Definition: ConvergenceReport.hpp:38
Definition: DeferredLogger.hpp:57
void debug(const std::string &tag, const std::string &message)
Definition: ISTLSolver.hpp:149
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:34
Vector extractVector(const Vector &x, const std::vector< int > &indices)
Definition: extractMatrix.hpp:104
void copySubMatrix(const Matrix &A, const std::vector< int > &indices, Matrix &B)
Definition: extractMatrix.hpp:35
void setGlobal(const Vector &x, const std::vector< int > &indices, Vector &global_x)
Definition: extractMatrix.hpp:115
Matrix extractMatrix(const Matrix &m, const std::vector< int > &indices)
Definition: extractMatrix.hpp:47
void unPack(PV &privar, const std::size_t meanings)
Definition: priVarsPacking.hpp:41
std::size_t pack(const PV &privar)
Definition: priVarsPacking.hpp:31
std::size_t countGlobalCells(const Grid &grid)
Get the number of cells of a global grid.
Definition: countGlobalCells.hpp:80
void detectOscillations(const std::vector< std::vector< Scalar > > &residualHistory, const int it, const int numPhases, const Scalar relaxRelTol, const int minimumOscillatingPhases, bool &oscillate, bool &stagnate)
Detect oscillation or stagnation in a given residual history.
Definition: blackoilboundaryratevector.hh:39
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Parallel::Communication communicator)
Create a global log combining local logs.
std::pair< std::vector< int >, int > partitionCells(const std::string &method, const int num_local_domains, const GridView &grid_view, const std::vector< Well > &wells, const std::unordered_map< std::string, std::set< int > > &possibleFutureConnections, const ZoltanPartitioningControl< Element > &zoltan_ctrl, const int num_neighbor_levels)
void printDomainDistributionSummary(const std::vector< int > &partition_vector, const std::vector< Domain > &domains, SimulatorReport &local_reports_accumulated, std::vector< SimulatorReport > &domain_reports_accumulated, const Grid &grid, int num_wells)
Definition: NlddReporting.hpp:219
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
void writePartitions(const std::filesystem::path &odir, const std::vector< Domain > &domains, const Grid &grid, const ElementMapper &elementMapper, const CartMapper &cartMapper)
Definition: NlddReporting.hpp:164
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
void writeNonlinearIterationsPerCell(const std::filesystem::path &odir, const std::vector< Domain > &domains, const std::vector< SimulatorReport > &domain_reports, const Grid &grid, const ElementMapper &elementMapper, const CartMapper &cartMapper)
Definition: NlddReporting.hpp:104
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:179
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition: FlowLinearSolverParameters.hpp:97
bool is_nldd_local_solver_
Definition: FlowLinearSolverParameters.hpp:111
void init(bool cprRequestedInDataFile)
bool linear_solver_print_json_definition_
Definition: FlowLinearSolverParameters.hpp:113
std::string linsolver_
Definition: FlowLinearSolverParameters.hpp:112
Definition: SimulatorReport.hpp:122
SimulatorReportSingle success
Definition: SimulatorReport.hpp:123
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34
double linear_solve_time
Definition: SimulatorReport.hpp:43
int skipped_domains
Definition: SimulatorReport.hpp:71
double assemble_time
Definition: SimulatorReport.hpp:39
double assemble_time_well
Definition: SimulatorReport.hpp:41
double solver_time
Definition: SimulatorReport.hpp:38
bool converged
Definition: SimulatorReport.hpp:55
double pre_post_time
Definition: SimulatorReport.hpp:40
double total_time
Definition: SimulatorReport.hpp:37
double linear_solve_setup_time
Definition: SimulatorReport.hpp:42
unsigned int total_newton_iterations
Definition: SimulatorReport.hpp:50
double update_time
Definition: SimulatorReport.hpp:45
double local_solve_time
Definition: SimulatorReport.hpp:44
unsigned int total_linearizations
Definition: SimulatorReport.hpp:49
unsigned int total_linear_iterations
Definition: SimulatorReport.hpp:51
Definition: SubDomain.hpp:85