My Project
Loading...
Searching...
No Matches
BlackoilModelNldd.hpp
1/*
2 Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3 Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4 Copyright 2014, 2015 Statoil ASA.
5 Copyright 2015 NTNU
6 Copyright 2015, 2016, 2017 IRIS AS
7
8 This file is part of the Open Porous Media project (OPM).
9
10 OPM is free software: you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14
15 OPM is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with OPM. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#ifndef OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED
26
27#include <dune/common/timer.hh>
28
29#include <opm/grid/common/SubGridPart.hpp>
30
31#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
32
33#include <opm/simulators/flow/countGlobalCells.hpp>
34#include <opm/simulators/flow/partitionCells.hpp>
35#include <opm/simulators/flow/priVarsPacking.hpp>
36#include <opm/simulators/flow/NonlinearSolver.hpp>
37#include <opm/simulators/flow/SubDomain.hpp>
38
39#include <opm/simulators/linalg/extractMatrix.hpp>
40
41#if COMPILE_BDA_BRIDGE
42#include <opm/simulators/linalg/ISTLSolverBda.hpp>
43#else
44#include <opm/simulators/linalg/ISTLSolver.hpp>
45#endif
46
47#include <opm/simulators/timestepping/ConvergenceReport.hpp>
48#include <opm/simulators/timestepping/SimulatorReport.hpp>
49#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
50
51#include <opm/simulators/utils/ComponentName.hpp>
52
53#include <fmt/format.h>
54
55#include <algorithm>
56#include <array>
57#include <cassert>
58#include <cmath>
59#include <cstddef>
60#include <filesystem>
61#include <fstream>
62#include <functional>
63#include <iomanip>
64#include <ios>
65#include <memory>
66#include <numeric>
67#include <sstream>
68#include <string>
69#include <type_traits>
70#include <utility>
71#include <vector>
72
73namespace Opm {
74
75template<class TypeTag> class BlackoilModel;
76
78template <class TypeTag>
80public:
88
89 using BVector = typename BlackoilModel<TypeTag>::BVector;
90 using Domain = SubDomain<Grid>;
92 using Mat = typename BlackoilModel<TypeTag>::Mat;
93
94 static constexpr int numEq = Indices::numEq;
95
101 : model_(model), rank_(model_.simulator().vanguard().grid().comm().rank())
102 {
103 // Create partitions.
104 const auto& [partition_vector, num_domains] = this->partitionCells();
105
106 // Scan through partitioning to get correct size for each.
107 std::vector<int> sizes(num_domains, 0);
108 for (const auto& p : partition_vector) {
109 ++sizes[p];
110 }
111
112 // Set up correctly sized vectors of entity seeds and of indices for each partition.
113 using EntitySeed = typename Grid::template Codim<0>::EntitySeed;
114 std::vector<std::vector<EntitySeed>> seeds(num_domains);
115 std::vector<std::vector<int>> partitions(num_domains);
116 for (int domain = 0; domain < num_domains; ++domain) {
117 seeds[domain].resize(sizes[domain]);
118 partitions[domain].resize(sizes[domain]);
119 }
120
121 // Iterate through grid once, setting the seeds of all partitions.
122 // Note: owned cells only!
123 const auto& grid = model_.simulator().vanguard().grid();
124
125 std::vector<int> count(num_domains, 0);
126 const auto& gridView = grid.leafGridView();
127 const auto beg = gridView.template begin<0, Dune::Interior_Partition>();
128 const auto end = gridView.template end<0, Dune::Interior_Partition>();
129 int cell = 0;
130 for (auto it = beg; it != end; ++it, ++cell) {
131 const int p = partition_vector[cell];
132 seeds[p][count[p]] = it->seed();
133 partitions[p][count[p]] = cell;
134 ++count[p];
135 }
136 assert(count == sizes);
137
138 // Create the domains.
139 for (int index = 0; index < num_domains; ++index) {
140 std::vector<bool> interior(partition_vector.size(), false);
141 for (int ix : partitions[index]) {
142 interior[ix] = true;
143 }
144
145 Dune::SubGridPart<Grid> view{grid, std::move(seeds[index])};
146
147 this->domains_.emplace_back(index,
148 std::move(partitions[index]),
149 std::move(interior),
150 std::move(view));
151 }
152
153 // Set up container for the local system matrices.
154 domain_matrices_.resize(num_domains);
155
156 // Set up container for the local linear solvers.
157 for (int index = 0; index < num_domains; ++index) {
158 // TODO: The ISTLSolver constructor will make
159 // parallel structures appropriate for the full grid
160 // only. This must be addressed before going parallel.
161 const auto& eclState = model_.simulator().vanguard().eclState();
163 loc_param.is_nldd_local_solver_ = true;
164 loc_param.init(eclState.getSimulationConfig().useCPR());
165 // Override solver type with umfpack if small domain.
166 if (domains_[index].cells.size() < 200) {
167 loc_param.linsolver_ = "umfpack";
168 }
169 loc_param.linear_solver_print_json_definition_ = false;
170 const bool force_serial = true;
171 domain_linsolvers_.emplace_back(model_.simulator(), loc_param, force_serial);
172 domain_linsolvers_.back().setDomainIndex(index);
173 }
174
175 assert(int(domains_.size()) == num_domains);
176 }
177
180 {
181 // Setup domain->well mapping.
182 model_.wellModel().setupDomains(domains_);
183 }
184
186 template <class NonlinearSolverType>
188 const SimulatorTimerInterface& timer,
189 NonlinearSolverType& nonlinear_solver)
190 {
191 // ----------- Set up reports and timer -----------
193 Dune::Timer perfTimer;
194
195 model_.initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
196
197 if (report.converged) {
198 return report;
199 }
200
201 // ----------- If not converged, do an NLDD iteration -----------
202
203 auto& solution = model_.simulator().model().solution(0);
204 auto initial_solution = solution;
205 auto locally_solved = initial_solution;
206
207 // ----------- Decide on an ordering for the domains -----------
208 const auto domain_order = this->getSubdomainOrder();
209
210 // ----------- Solve each domain separately -----------
211 DeferredLogger logger;
212 std::vector<SimulatorReportSingle> domain_reports(domains_.size());
213 for (const int domain_index : domain_order) {
214 const auto& domain = domains_[domain_index];
215 SimulatorReportSingle local_report;
216 try {
217 switch (model_.param().local_solve_approach_) {
218 case DomainSolveApproach::Jacobi:
219 solveDomainJacobi(solution, locally_solved, local_report, logger,
220 iteration, timer, domain);
221 break;
222 default:
223 case DomainSolveApproach::GaussSeidel:
224 solveDomainGaussSeidel(solution, locally_solved, local_report, logger,
225 iteration, timer, domain);
226 break;
227 }
228 }
229 catch (...) {
230 // Something went wrong during local solves.
231 local_report.converged = false;
232 }
233 // This should have updated the global matrix to be
234 // dR_i/du_j evaluated at new local solutions for
235 // i == j, at old solution for i != j.
236 if (!local_report.converged) {
237 // TODO: more proper treatment, including in parallel.
238 logger.debug(fmt::format("Convergence failure in domain {} on rank {}." , domain.index, rank_));
239 }
240 domain_reports[domain.index] = local_report;
241 }
242
243 // Communicate and log all messages.
244 auto global_logger = gatherDeferredLogger(logger, model_.simulator().vanguard().grid().comm());
245 global_logger.logMessages();
246
247 // Accumulate local solve data.
248 // Putting the counts in a single array to avoid multiple
249 // comm.sum() calls. Keeping the named vars for readability.
250 std::array<int, 4> counts{ 0, 0, 0, static_cast<int>(domain_reports.size()) };
251 int& num_converged = counts[0];
252 int& num_converged_already = counts[1];
253 int& num_local_newtons = counts[2];
254 int& num_domains = counts[3];
255 {
257 for (const auto& dr : domain_reports) {
258 if (dr.converged) {
259 ++num_converged;
260 if (dr.total_newton_iterations == 0) {
261 ++num_converged_already;
262 }
263 }
264 rep += dr;
265 }
266 num_local_newtons = rep.total_newton_iterations;
267 local_reports_accumulated_ += rep;
268 }
269
270 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
271 solution = locally_solved;
272 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
273 }
274
275#if HAVE_MPI
276 // Communicate solutions:
277 // With multiple processes, this process' overlap (i.e. not
278 // owned) cells' solution values have been modified by local
279 // solves in the owning processes, and remain unchanged
280 // here. We must therefore receive the updated solution on the
281 // overlap cells and update their intensive quantities before
282 // we move on.
283 const auto& comm = model_.simulator().vanguard().grid().comm();
284 if (comm.size() > 1) {
285 const auto* ccomm = model_.simulator().model().newtonMethod().linearSolver().comm();
286
287 // Copy numerical values from primary vars.
288 ccomm->copyOwnerToAll(solution, solution);
289
290 // Copy flags from primary vars.
291 const std::size_t num = solution.size();
292 Dune::BlockVector<std::size_t> allmeanings(num);
293 for (std::size_t ii = 0; ii < num; ++ii) {
294 allmeanings[ii] = PVUtil::pack(solution[ii]);
295 }
296 ccomm->copyOwnerToAll(allmeanings, allmeanings);
297 for (std::size_t ii = 0; ii < num; ++ii) {
298 PVUtil::unPack(solution[ii], allmeanings[ii]);
299 }
300
301 // Update intensive quantities for our overlap values.
302 model_.simulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(/*timeIdx=*/0);
303
304 // Make total counts of domains converged.
305 comm.sum(counts.data(), counts.size());
306 }
307#endif // HAVE_MPI
308
309 const bool is_iorank = this->rank_ == 0;
310 if (is_iorank) {
311 OpmLog::debug(fmt::format("Local solves finished. Converged for {}/{} domains. {} domains did no work. {} total local Newton iterations.\n",
312 num_converged, num_domains, num_converged_already, num_local_newtons));
313 }
314
315 // Finish with a Newton step.
316 // Note that the "iteration + 100" is a simple way to avoid entering
317 // "if (iteration == 0)" and similar blocks, and also makes it a little
318 // easier to spot the iteration residuals in the DBG file. A more sophisticated
319 // approach can be done later.
320 auto rep = model_.nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver);
321 report += rep;
322 if (rep.converged) {
323 report.converged = true;
324 }
325 return report;
326 }
327
330 {
331 return local_reports_accumulated_;
332 }
333
334 void writePartitions(const std::filesystem::path& odir) const
335 {
336 const auto& elementMapper = this->model_.simulator().model().elementMapper();
337 const auto& cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
338
339 const auto& grid = this->model_.simulator().vanguard().grid();
340 const auto& comm = grid.comm();
341 const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
342
343 std::ofstream pfile { odir / fmt::format("{1:0>{0}}", nDigit, comm.rank()) };
344
345 const auto p = this->reconstitutePartitionVector();
346 auto i = 0;
347 for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
348 pfile << comm.rank() << ' '
349 << cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
350 << p[i++] << '\n';
351 }
352 }
353
354private:
355
357 std::pair<SimulatorReportSingle, ConvergenceReport>
358 solveDomain(const Domain& domain,
359 const SimulatorTimerInterface& timer,
360 DeferredLogger& logger,
361 [[maybe_unused]] const int global_iteration,
362 const bool initial_assembly_required)
363 {
364 auto& modelSimulator = model_.simulator();
365
366 SimulatorReportSingle report;
367 Dune::Timer solveTimer;
368 solveTimer.start();
369 Dune::Timer detailTimer;
370
371 modelSimulator.model().newtonMethod().setIterationIndex(0);
372
373 // When called, if assembly has already been performed
374 // with the initial values, we only need to check
375 // for local convergence. Otherwise, we must do a local
376 // assembly.
377 int iter = 0;
378 if (initial_assembly_required) {
379 detailTimer.start();
380 modelSimulator.model().newtonMethod().setIterationIndex(iter);
381 // TODO: we should have a beginIterationLocal function()
382 // only handling the well model for now
383 modelSimulator.problem().wellModel().assembleDomain(modelSimulator.model().newtonMethod().numIterations(),
384 modelSimulator.timeStepSize(),
385 domain);
386 // Assemble reservoir locally.
387 report += this->assembleReservoirDomain(domain);
388 report.assemble_time += detailTimer.stop();
389 }
390 detailTimer.reset();
391 detailTimer.start();
392 std::vector<Scalar> resnorms;
393 auto convreport = this->getDomainConvergence(domain, timer, 0, logger, resnorms);
394 if (convreport.converged()) {
395 // TODO: set more info, timing etc.
396 report.converged = true;
397 return { report, convreport };
398 }
399
400 // We have already assembled for the first iteration,
401 // but not done the Schur complement for the wells yet.
402 detailTimer.reset();
403 detailTimer.start();
404 model_.wellModel().linearizeDomain(domain,
405 modelSimulator.model().linearizer().jacobian(),
406 modelSimulator.model().linearizer().residual());
407 const double tt1 = detailTimer.stop();
408 report.assemble_time += tt1;
409 report.assemble_time_well += tt1;
410
411 // Local Newton loop.
412 const int max_iter = model_.param().max_local_solve_iterations_;
413 const auto& grid = modelSimulator.vanguard().grid();
414 double damping_factor = 1.0;
415 std::vector<std::vector<Scalar>> convergence_history;
416 convergence_history.reserve(20);
417 convergence_history.push_back(resnorms);
418 do {
419 // Solve local linear system.
420 // Note that x has full size, we expect it to be nonzero only for in-domain cells.
421 const int nc = grid.size(0);
422 BVector x(nc);
423 detailTimer.reset();
424 detailTimer.start();
425 this->solveJacobianSystemDomain(domain, x);
426 model_.wellModel().postSolveDomain(x, domain);
427 if (damping_factor != 1.0) {
428 x *= damping_factor;
429 }
430 report.linear_solve_time += detailTimer.stop();
431 report.linear_solve_setup_time += model_.linearSolveSetupTime();
432 report.total_linear_iterations = model_.linearIterationsLastSolve();
433
434 // Update local solution. // TODO: x is still full size, should we optimize it?
435 detailTimer.reset();
436 detailTimer.start();
437 this->updateDomainSolution(domain, x);
438 report.update_time += detailTimer.stop();
439
440 // Assemble well and reservoir.
441 detailTimer.reset();
442 detailTimer.start();
443 ++iter;
444 modelSimulator.model().newtonMethod().setIterationIndex(iter);
445 // TODO: we should have a beginIterationLocal function()
446 // only handling the well model for now
447 // Assemble reservoir locally.
448 modelSimulator.problem().wellModel().assembleDomain(modelSimulator.model().newtonMethod().numIterations(),
449 modelSimulator.timeStepSize(),
450 domain);
451 report += this->assembleReservoirDomain(domain);
452 report.assemble_time += detailTimer.stop();
453
454 // Check for local convergence.
455 detailTimer.reset();
456 detailTimer.start();
457 resnorms.clear();
458 convreport = this->getDomainConvergence(domain, timer, iter, logger, resnorms);
459 convergence_history.push_back(resnorms);
460
461 // apply the Schur complement of the well model to the
462 // reservoir linearized equations
463 detailTimer.reset();
464 detailTimer.start();
465 model_.wellModel().linearizeDomain(domain,
466 modelSimulator.model().linearizer().jacobian(),
467 modelSimulator.model().linearizer().residual());
468 const double tt2 = detailTimer.stop();
469 report.assemble_time += tt2;
470 report.assemble_time_well += tt2;
471
472 // Check if we should dampen. Only do so if wells are converged.
473 if (!convreport.converged() && !convreport.wellFailed()) {
474 bool oscillate = false;
475 bool stagnate = false;
476 const int numPhases = convergence_history.front().size();
477 detail::detectOscillations(convergence_history, iter, numPhases,
478 Scalar{0.2}, 1, oscillate, stagnate);
479 if (oscillate) {
480 damping_factor *= 0.85;
481 logger.debug(fmt::format("| Damping factor is now {}", damping_factor));
482 }
483 }
484 } while (!convreport.converged() && iter <= max_iter);
485
486 modelSimulator.problem().endIteration();
487
488 report.converged = convreport.converged();
489 report.total_newton_iterations = iter;
490 report.total_linearizations = iter;
491 report.total_time = solveTimer.stop();
492 // TODO: set more info, timing etc.
493 return { report, convreport };
494 }
495
497 SimulatorReportSingle assembleReservoirDomain(const Domain& domain)
498 {
499 // -------- Mass balance equations --------
500 model_.simulator().model().linearizer().linearizeDomain(domain);
501 return model_.wellModel().lastReport();
502 }
503
505 void solveJacobianSystemDomain(const Domain& domain, BVector& global_x)
506 {
507 const auto& modelSimulator = model_.simulator();
508
509 Dune::Timer perfTimer;
510 perfTimer.start();
511
512 const Mat& main_matrix = modelSimulator.model().linearizer().jacobian().istlMatrix();
513 if (domain_matrices_[domain.index]) {
514 Details::copySubMatrix(main_matrix, domain.cells, *domain_matrices_[domain.index]);
515 } else {
516 domain_matrices_[domain.index] = std::make_unique<Mat>(Details::extractMatrix(main_matrix, domain.cells));
517 }
518 auto& jac = *domain_matrices_[domain.index];
519 auto res = Details::extractVector(modelSimulator.model().linearizer().residual(),
520 domain.cells);
521 auto x = res;
522
523 // set initial guess
524 global_x = 0.0;
525 x = 0.0;
526
527 auto& linsolver = domain_linsolvers_[domain.index];
528
529 linsolver.prepare(jac, res);
530 model_.linearSolveSetupTime() = perfTimer.stop();
531 linsolver.setResidual(res);
532 linsolver.solve(x);
533
534 Details::setGlobal(x, domain.cells, global_x);
535 }
536
538 void updateDomainSolution(const Domain& domain, const BVector& dx)
539 {
540 auto& simulator = model_.simulator();
541 auto& newtonMethod = simulator.model().newtonMethod();
542 SolutionVector& solution = simulator.model().solution(/*timeIdx=*/0);
543
544 newtonMethod.update_(/*nextSolution=*/solution,
545 /*curSolution=*/solution,
546 /*update=*/dx,
547 /*resid=*/dx,
548 domain.cells); // the update routines of the black
549 // oil model do not care about the
550 // residual
551
552 // if the solution is updated, the intensive quantities need to be recalculated
553 simulator.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
554 }
555
557 std::pair<Scalar, Scalar> localDomainConvergenceData(const Domain& domain,
558 std::vector<Scalar>& R_sum,
559 std::vector<Scalar>& maxCoeff,
560 std::vector<Scalar>& B_avg,
561 std::vector<int>& maxCoeffCell)
562 {
563 const auto& modelSimulator = model_.simulator();
564
565 Scalar pvSumLocal = 0.0;
566 Scalar numAquiferPvSumLocal = 0.0;
567 const auto& model = modelSimulator.model();
568 const auto& problem = modelSimulator.problem();
569
570 const auto& modelResid = modelSimulator.model().linearizer().residual();
571
572 ElementContext elemCtx(modelSimulator);
573 const auto& gridView = domain.view;
574 const auto& elemEndIt = gridView.template end</*codim=*/0>();
575 IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
576
577 for (auto elemIt = gridView.template begin</*codim=*/0>();
578 elemIt != elemEndIt;
579 ++elemIt)
580 {
581 if (elemIt->partitionType() != Dune::InteriorEntity) {
582 continue;
583 }
584 const auto& elem = *elemIt;
585 elemCtx.updatePrimaryStencil(elem);
586 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
587
588 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
589 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
590 const auto& fs = intQuants.fluidState();
591
592 const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
593 model.dofTotalVolume(cell_idx);
594 pvSumLocal += pvValue;
595
596 if (isNumericalAquiferCell(elem))
597 {
598 numAquiferPvSumLocal += pvValue;
599 }
600
601 model_.getMaxCoeff(cell_idx, intQuants, fs, modelResid, pvValue,
602 B_avg, R_sum, maxCoeff, maxCoeffCell);
603 }
604
605 // compute local average in terms of global number of elements
606 const int bSize = B_avg.size();
607 for ( int i = 0; i<bSize; ++i )
608 {
609 B_avg[ i ] /= Scalar(domain.cells.size());
610 }
611
612 return {pvSumLocal, numAquiferPvSumLocal};
613 }
614
615 ConvergenceReport getDomainReservoirConvergence(const double reportTime,
616 const double dt,
617 const int iteration,
618 const Domain& domain,
619 DeferredLogger& logger,
620 std::vector<Scalar>& B_avg,
621 std::vector<Scalar>& residual_norms)
622 {
623 using Vector = std::vector<Scalar>;
624
625 const int numComp = numEq;
626 Vector R_sum(numComp, 0.0 );
627 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest() );
628 std::vector<int> maxCoeffCell(numComp, -1);
629 const auto [ pvSum, numAquiferPvSum]
630 = this->localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell);
631
632 auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt);
633 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
634
635 // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
636 // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
637 // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
638 const bool use_relaxed_cnv = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
639 iteration >= model_.param().min_strict_cnv_iter_;
640 // Tighter bound for local convergence should increase the
641 // likelyhood of: local convergence => global convergence
642 const Scalar tol_cnv = model_.param().local_tolerance_scaling_cnv_
643 * (use_relaxed_cnv ? model_.param().tolerance_cnv_relaxed_
644 : model_.param().tolerance_cnv_);
645
646 const bool use_relaxed_mb = iteration >= model_.param().min_strict_mb_iter_;
647 const Scalar tol_mb = model_.param().local_tolerance_scaling_mb_
648 * (use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
649
650 // Finish computation
651 std::vector<Scalar> CNV(numComp);
652 std::vector<Scalar> mass_balance_residual(numComp);
653 for (int compIdx = 0; compIdx < numComp; ++compIdx )
654 {
655 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
656 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
657 residual_norms.push_back(CNV[compIdx]);
658 }
659
660 // Create convergence report.
661 ConvergenceReport report{reportTime};
662 using CR = ConvergenceReport;
663 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
664 Scalar res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
665 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
666 CR::ReservoirFailure::Type::Cnv };
667 Scalar tol[2] = { tol_mb, tol_cnv };
668 for (int ii : {0, 1}) {
669 if (std::isnan(res[ii])) {
670 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
671 logger.debug("NaN residual for " + model_.compNames().name(compIdx) + " equation.");
672 } else if (res[ii] > model_.param().max_residual_allowed_) {
673 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
674 logger.debug("Too large residual for " + model_.compNames().name(compIdx) + " equation.");
675 } else if (res[ii] < 0.0) {
676 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
677 logger.debug("Negative residual for " + model_.compNames().name(compIdx) + " equation.");
678 } else if (res[ii] > tol[ii]) {
679 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
680 }
681
682 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
683 }
684 }
685
686 // Output of residuals. If converged at initial state, log nothing.
687 const bool converged_at_initial_state = (report.converged() && iteration == 0);
688 if (!converged_at_initial_state) {
689 if (iteration == 0) {
690 // Log header.
691 std::string msg = fmt::format("Domain {} on rank {}, size {}, containing cell {}\n| Iter",
692 domain.index, this->rank_, domain.cells.size(), domain.cells[0]);
693 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
694 msg += " MB(";
695 msg += model_.compNames().name(compIdx)[0];
696 msg += ") ";
697 }
698 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
699 msg += " CNV(";
700 msg += model_.compNames().name(compIdx)[0];
701 msg += ") ";
702 }
703 logger.debug(msg);
704 }
705 // Log convergence data.
706 std::ostringstream ss;
707 ss << "| ";
708 const std::streamsize oprec = ss.precision(3);
709 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
710 ss << std::setw(4) << iteration;
711 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
712 ss << std::setw(11) << mass_balance_residual[compIdx];
713 }
714 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
715 ss << std::setw(11) << CNV[compIdx];
716 }
717 ss.precision(oprec);
718 ss.flags(oflags);
719 logger.debug(ss.str());
720 }
721
722 return report;
723 }
724
725 ConvergenceReport getDomainConvergence(const Domain& domain,
726 const SimulatorTimerInterface& timer,
727 const int iteration,
728 DeferredLogger& logger,
729 std::vector<Scalar>& residual_norms)
730 {
731 std::vector<Scalar> B_avg(numEq, 0.0);
732 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
733 timer.currentStepLength(),
734 iteration,
735 domain,
736 logger,
737 B_avg,
738 residual_norms);
739 report += model_.wellModel().getDomainWellConvergence(domain, B_avg, logger);
740 return report;
741 }
742
744 std::vector<int> getSubdomainOrder()
745 {
746 const auto& modelSimulator = model_.simulator();
747 const auto& solution = modelSimulator.model().solution(0);
748
749 std::vector<int> domain_order(domains_.size());
750 std::iota(domain_order.begin(), domain_order.end(), 0);
751
752 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
753 // Do nothing, 0..n-1 order is fine.
754 return domain_order;
755 } else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) {
756 // Calculate the measure used to order the domains.
757 std::vector<Scalar> measure_per_domain(domains_.size());
758 switch (model_.param().local_domain_ordering_) {
759 case DomainOrderingMeasure::AveragePressure: {
760 // Use average pressures to order domains.
761 for (const auto& domain : domains_) {
762 Scalar press_sum = 0.0;
763 for (const int c : domain.cells) {
764 press_sum += solution[c][Indices::pressureSwitchIdx];
765 }
766 const Scalar avgpress = press_sum / domain.cells.size();
767 measure_per_domain[domain.index] = avgpress;
768 }
769 break;
770 }
771 case DomainOrderingMeasure::MaxPressure: {
772 // Use max pressures to order domains.
773 for (const auto& domain : domains_) {
774 Scalar maxpress = 0.0;
775 for (const int c : domain.cells) {
776 maxpress = std::max(maxpress, solution[c][Indices::pressureSwitchIdx]);
777 }
778 measure_per_domain[domain.index] = maxpress;
779 }
780 break;
781 }
782 case DomainOrderingMeasure::Residual: {
783 // Use maximum residual to order domains.
784 const auto& residual = modelSimulator.model().linearizer().residual();
785 const int num_vars = residual[0].size();
786 for (const auto& domain : domains_) {
787 Scalar maxres = 0.0;
788 for (const int c : domain.cells) {
789 for (int ii = 0; ii < num_vars; ++ii) {
790 maxres = std::max(maxres, std::fabs(residual[c][ii]));
791 }
792 }
793 measure_per_domain[domain.index] = maxres;
794 }
795 break;
796 }
797 } // end of switch (model_.param().local_domain_ordering_)
798
799 // Sort by largest measure, keeping index order if equal.
800 const auto& m = measure_per_domain;
801 std::stable_sort(domain_order.begin(), domain_order.end(),
802 [&m](const int i1, const int i2){ return m[i1] > m[i2]; });
803 return domain_order;
804 } else {
805 throw std::logic_error("Domain solve approach must be Jacobi or Gauss-Seidel");
806 }
807 }
808
809 template<class GlobalEqVector>
810 void solveDomainJacobi(GlobalEqVector& solution,
811 GlobalEqVector& locally_solved,
812 SimulatorReportSingle& local_report,
813 DeferredLogger& logger,
814 const int iteration,
815 const SimulatorTimerInterface& timer,
816 const Domain& domain)
817 {
818 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
819 auto initial_local_solution = Details::extractVector(solution, domain.cells);
820 auto res = solveDomain(domain, timer, logger, iteration, false);
821 local_report = res.first;
822 if (local_report.converged) {
823 auto local_solution = Details::extractVector(solution, domain.cells);
824 Details::setGlobal(local_solution, domain.cells, locally_solved);
825 Details::setGlobal(initial_local_solution, domain.cells, solution);
826 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
827 } else {
828 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
829 Details::setGlobal(initial_local_solution, domain.cells, solution);
830 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
831 }
832 }
833
834 template<class GlobalEqVector>
835 void solveDomainGaussSeidel(GlobalEqVector& solution,
836 GlobalEqVector& locally_solved,
837 SimulatorReportSingle& local_report,
838 DeferredLogger& logger,
839 const int iteration,
840 const SimulatorTimerInterface& timer,
841 const Domain& domain)
842 {
843 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
844 auto initial_local_solution = Details::extractVector(solution, domain.cells);
845 auto res = solveDomain(domain, timer, logger, iteration, true);
846 local_report = res.first;
847 if (!local_report.converged) {
848 // We look at the detailed convergence report to evaluate
849 // if we should accept the unconverged solution.
850 const auto& convrep = res.second;
851 // We do not accept a solution if the wells are unconverged.
852 if (!convrep.wellFailed()) {
853 // Calculare the sums of the mb and cnv failures.
854 Scalar mb_sum = 0.0;
855 Scalar cnv_sum = 0.0;
856 for (const auto& rc : convrep.reservoirConvergence()) {
857 if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
858 mb_sum += rc.value();
859 } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) {
860 cnv_sum += rc.value();
861 }
862 }
863 // If not too high, we overrule the convergence failure.
864 const Scalar acceptable_local_mb_sum = 1e-3;
865 const Scalar acceptable_local_cnv_sum = 1.0;
866 if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) {
867 local_report.converged = true;
868 logger.debug(fmt::format("Accepting solution in unconverged domain {} on rank {}.", domain.index, rank_));
869 logger.debug(fmt::format("Value of mb_sum: {} cnv_sum: {}", mb_sum, cnv_sum));
870 } else {
871 logger.debug("Unconverged local solution.");
872 }
873 } else {
874 logger.debug("Unconverged local solution with well convergence failures:");
875 for (const auto& wf : convrep.wellFailures()) {
876 logger.debug(to_string(wf));
877 }
878 }
879 }
880 if (local_report.converged) {
881 auto local_solution = Details::extractVector(solution, domain.cells);
882 Details::setGlobal(local_solution, domain.cells, locally_solved);
883 } else {
884 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
885 Details::setGlobal(initial_local_solution, domain.cells, solution);
886 model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
887 }
888 }
889
890 Scalar computeCnvErrorPvLocal(const Domain& domain,
891 const std::vector<Scalar>& B_avg, double dt) const
892 {
893 Scalar errorPV{};
894 const auto& simulator = model_.simulator();
895 const auto& model = simulator.model();
896 const auto& problem = simulator.problem();
897 const auto& residual = simulator.model().linearizer().residual();
898
899 for (const int cell_idx : domain.cells) {
900 const Scalar pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
901 model.dofTotalVolume(cell_idx);
902 const auto& cellResidual = residual[cell_idx];
903 bool cnvViolated = false;
904
905 for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) {
906 using std::fabs;
907 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
908 cnvViolated = cnvViolated || (fabs(CNV) > model_.param().tolerance_cnv_);
909 }
910
911 if (cnvViolated) {
912 errorPV += pvValue;
913 }
914 }
915 return errorPV;
916 }
917
918 decltype(auto) partitionCells() const
919 {
920 const auto& grid = this->model_.simulator().vanguard().grid();
921
922 using GridView = std::remove_cv_t<std::remove_reference_t<
923 decltype(grid.leafGridView())>>;
924
925 using Element = std::remove_cv_t<std::remove_reference_t<
926 typename GridView::template Codim<0>::Entity>>;
927
928 const auto& param = this->model_.param();
929
930 auto zoltan_ctrl = ZoltanPartitioningControl<Element>{};
931
932 zoltan_ctrl.domain_imbalance = param.local_domain_partition_imbalance_;
933
934 zoltan_ctrl.index =
935 [elementMapper = &this->model_.simulator().model().elementMapper()]
936 (const Element& element)
937 {
938 return elementMapper->index(element);
939 };
940
941 zoltan_ctrl.local_to_global =
942 [cartMapper = &this->model_.simulator().vanguard().cartesianIndexMapper()]
943 (const int elemIdx)
944 {
945 return cartMapper->cartesianIndex(elemIdx);
946 };
947
948 // Forming the list of wells is expensive, so do this only if needed.
949 const auto need_wells = param.local_domain_partition_method_ == "zoltan";
950
951 const auto wells = need_wells
952 ? this->model_.simulator().vanguard().schedule().getWellsatEnd()
953 : std::vector<Well>{};
954
955 const auto& possibleFutureConnectionSet = need_wells
956 ? this->model_.simulator().vanguard().schedule().getPossibleFutureConnections()
957 : std::unordered_map<std::string, std::set<int>> {};
958
959 // If defaulted parameter for number of domains, choose a reasonable default.
960 constexpr int default_cells_per_domain = 1000;
961 const int num_domains = (param.num_local_domains_ > 0)
962 ? param.num_local_domains_
963 : detail::countGlobalCells(grid) / default_cells_per_domain;
964
965 return ::Opm::partitionCells(param.local_domain_partition_method_,
966 num_domains, grid.leafGridView(), wells,
967 possibleFutureConnectionSet, zoltan_ctrl);
968 }
969
970 std::vector<int> reconstitutePartitionVector() const
971 {
972 const auto& grid = this->model_.simulator().vanguard().grid();
973
974 auto numD = std::vector<int>(grid.comm().size() + 1, 0);
975 numD[grid.comm().rank() + 1] = static_cast<int>(this->domains_.size());
976 grid.comm().sum(numD.data(), numD.size());
977 std::partial_sum(numD.begin(), numD.end(), numD.begin());
978
979 auto p = std::vector<int>(grid.size(0));
980 auto maxCellIdx = std::numeric_limits<int>::min();
981
982 auto d = numD[grid.comm().rank()];
983 for (const auto& domain : this->domains_) {
984 for (const auto& cell : domain.cells) {
985 p[cell] = d;
986 if (cell > maxCellIdx) {
987 maxCellIdx = cell;
988 }
989 }
990
991 ++d;
992 }
993
994 p.erase(p.begin() + maxCellIdx + 1, p.end());
995 return p;
996 }
997
998 BlackoilModel<TypeTag>& model_;
999 std::vector<Domain> domains_;
1000 std::vector<std::unique_ptr<Mat>> domain_matrices_;
1001 std::vector<ISTLSolverType> domain_linsolvers_;
1002 SimulatorReportSingle local_reports_accumulated_;
1003 int rank_ = 0;
1004};
1005
1006} // namespace Opm
1007
1008#endif // OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED
A NLDD implementation for three-phase black oil.
Definition BlackoilModelNldd.hpp:79
BlackoilModelNldd(BlackoilModel< TypeTag > &model)
The constructor sets up the subdomains.
Definition BlackoilModelNldd.hpp:100
SimulatorReportSingle nonlinearIterationNldd(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Do one non-linear NLDD iteration.
Definition BlackoilModelNldd.hpp:187
const SimulatorReportSingle & localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModelNldd.hpp:329
void prepareStep()
Called before starting a time step.
Definition BlackoilModelNldd.hpp:179
A model implementation for three-phase black oil.
Definition BlackoilModelNldd.hpp:75
Definition DeferredLogger.hpp:57
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolver.hpp:144
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
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:235
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Opm::Parallel::Communication)
Create a global log combining local logs.
Definition gatherDeferredLogger.cpp:168
Solver parameters for the BlackoilModel.
Definition BlackoilModelParameters.hpp:162
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:95
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34
Representing a part of a grid, in a way suitable for performing local solves.
Definition SubDomain.hpp:62