My Project
Loading...
Searching...
No Matches
BlackoilModelEbosNldd.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_BLACKOILMODELEBOS_NLDD_HEADER_INCLUDED
25#define OPM_BLACKOILMODELEBOS_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/partitionCells.hpp>
34#include <opm/simulators/flow/priVarsPacking.hpp>
35#include <opm/simulators/flow/SubDomain.hpp>
36
37#include <opm/simulators/linalg/extractMatrix.hpp>
38
39#if COMPILE_BDA_BRIDGE
40#include <opm/simulators/linalg/ISTLSolverEbosBda.hpp>
41#else
42#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
43#endif
44
45#include <opm/simulators/timestepping/ConvergenceReport.hpp>
46#include <opm/simulators/timestepping/SimulatorReport.hpp>
47#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
48
49#include <opm/simulators/utils/ComponentName.hpp>
50
51#include <fmt/format.h>
52
53#include <algorithm>
54#include <cassert>
55#include <cmath>
56#include <cstddef>
57#include <iomanip>
58#include <ios>
59#include <memory>
60#include <numeric>
61#include <sstream>
62#include <string>
63#include <utility>
64#include <vector>
65
66namespace Opm {
67
68template<class TypeTag> class BlackoilModelEbos;
69
71template <class TypeTag>
73public:
81
82 using BVector = typename BlackoilModelEbos<TypeTag>::BVector;
83 using Domain = SubDomain<Grid>;
85 using Mat = typename BlackoilModelEbos<TypeTag>::Mat;
86
87 static constexpr int numEq = Indices::numEq;
88
94 : model_(model)
95 {
96 const auto& grid = model_.ebosSimulator().vanguard().grid();
97 const auto& schedule = model_.ebosSimulator().vanguard().schedule();
98
99 // Create partitions.
100 const auto& [partition_vector, num_domains] =
101 partitionCells(grid,
102 schedule.getWellsatEnd(),
103 model_.param().local_domain_partition_method_,
104 model_.param().num_local_domains_,
105 model_.param().local_domain_partition_imbalance_);
106
107 // Scan through partitioning to get correct size for each.
108 std::vector<int> sizes(num_domains, 0);
109 for (const auto& p : partition_vector) {
110 ++sizes[p];
111 }
112
113 // Set up correctly sized vectors of entity seeds and of indices for each partition.
114 using EntitySeed = typename Grid::template Codim<0>::EntitySeed;
115 std::vector<std::vector<EntitySeed>> seeds(num_domains);
116 std::vector<std::vector<int>> partitions(num_domains);
117 for (int domain = 0; domain < num_domains; ++domain) {
118 seeds[domain].resize(sizes[domain]);
119 partitions[domain].resize(sizes[domain]);
120 }
121
122 // Iterate through grid once, setting the seeds of all partitions.
123 // Note: owned cells only!
124 std::vector<int> count(num_domains, 0);
125 const auto& gridView = grid.leafGridView();
126 const auto beg = gridView.template begin<0, Dune::Interior_Partition>();
127 const auto end = gridView.template end<0, Dune::Interior_Partition>();
128 int cell = 0;
129 for (auto it = beg; it != end; ++it, ++cell) {
130 const int p = partition_vector[cell];
131 seeds[p][count[p]] = it->seed();
132 partitions[p][count[p]] = cell;
133 ++count[p];
134 }
135 assert(count == sizes);
136
137 // Create the domains.
138 for (int index = 0; index < num_domains; ++index) {
139 std::vector<bool> interior(partition_vector.size(), false);
140 for (int ix : partitions[index]) {
141 interior[ix] = true;
142 }
143
144 Dune::SubGridPart<Grid> view{grid, std::move(seeds[index])};
145
146 this->domains_.emplace_back(index,
147 std::move(partitions[index]),
148 std::move(interior),
149 std::move(view));
150 }
151
152 // Set up container for the local system matrices.
153 domain_matrices_.resize(num_domains);
154
155 // Set up container for the local linear solvers.
156 for (int index = 0; index < num_domains; ++index) {
157 // TODO: The ISTLSolverEbos constructor will make
158 // parallel structures appropriate for the full grid
159 // only. This must be addressed before going parallel.
160 const auto& eclState = model_.ebosSimulator().vanguard().eclState();
162 loc_param.template init<TypeTag>(eclState.getSimulationConfig().useCPR());
163 // Override solver type with umfpack if small domain.
164 // Otherwise hardcode to ILU0
165 if (domains_[index].cells.size() < 200) {
166 loc_param.linsolver_ = "umfpack";
167 } else {
168 loc_param.linsolver_ = "ilu0";
169 loc_param.linear_solver_reduction_ = 1e-2;
170 }
171 loc_param.linear_solver_print_json_definition_ = false;
172 const bool force_serial = true;
173 domain_linsolvers_.emplace_back(model_.ebosSimulator(), loc_param, force_serial);
174 }
175
176 assert(int(domains_.size()) == num_domains);
177 }
178
181 {
182 // Setup domain->well mapping.
183 model_.wellModel().setupDomains(domains_);
184 }
185
187 template <class NonlinearSolverType>
191 {
192 // ----------- Set up reports and timer -----------
194 Dune::Timer perfTimer;
195
196 model_.initialLinearization(report, iteration, nonlinear_solver.minIter(), timer);
197
198 if (report.converged) {
199 return report;
200 }
201
202 // ----------- If not converged, do an NLDD iteration -----------
203
204 auto& solution = model_.ebosSimulator().model().solution(0);
207
208 // ----------- Decide on an ordering for the domains -----------
209 const auto domain_order = this->getSubdomainOrder();
210
211 // ----------- Solve each domain separately -----------
212 std::vector<SimulatorReportSingle> domain_reports(domains_.size());
213 for (const int domain_index : domain_order) {
214 const auto& domain = domains_[domain_index];
216 switch (model_.param().local_solve_approach_) {
217 case DomainSolveApproach::Jacobi:
218 solveDomainJacobi(solution, locally_solved, local_report,
220 break;
221 default:
222 case DomainSolveApproach::GaussSeidel:
223 solveDomainGaussSeidel(solution, locally_solved, local_report,
225 break;
226 }
227 // This should have updated the global matrix to be
228 // dR_i/du_j evaluated at new local solutions for
229 // i == j, at old solution for i != j.
230 if (!local_report.converged) {
231 // TODO: more proper treatment, including in parallel.
232 OpmLog::debug("Convergence failure in domain " + std::to_string(domain.index));
233 }
235 }
236
237 // Log summary of local solve convergence to DBG file.
238 {
239 int num_converged = 0;
241 for (const auto& dr : domain_reports) {
242 if (dr.converged) {
244 }
245 rep += dr;
246 }
247 std::ostringstream os;
248 os << fmt::format("Local solves finished. Converged for {}/{} domains.\n",
250 rep.reportFullyImplicit(os, nullptr);
251 OpmLog::debug(os.str());
252 local_reports_accumulated_ += rep;
253 }
254
255 if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) {
257 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
258 }
259
260#if HAVE_MPI
261 // Communicate solutions:
262 // With multiple processes, this process' overlap (i.e. not
263 // owned) cells' solution values have been modified by local
264 // solves in the owning processes, and remain unchanged
265 // here. We must therefore receive the updated solution on the
266 // overlap cells and update their intensive quantities before
267 // we move on.
268 const auto& comm = model_.ebosSimulator().vanguard().grid().comm();
269 if (comm.size() > 1) {
270 const auto* ccomm = model_.ebosSimulator().model().newtonMethod().linearSolver().comm();
271
272 // Copy numerical values from primary vars.
273 ccomm->copyOwnerToAll(solution, solution);
274
275 // Copy flags from primary vars.
276 const std::size_t num = solution.size();
277 Dune::BlockVector<std::size_t> allmeanings(num);
278 for (std::size_t ii = 0; ii < num; ++ii) {
279 allmeanings[ii] = PVUtil::pack(solution[ii]);
280 }
281 ccomm->copyOwnerToAll(allmeanings, allmeanings);
282 for (std::size_t ii = 0; ii < num; ++ii) {
283 PVUtil::unPack(solution[ii], allmeanings[ii]);
284 }
285
286 // Update intensive quantities for our overlap values.
287 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(/*timeIdx=*/0);
288 }
289#endif // HAVE_MPI
290
291 // Finish with a Newton step.
292 // Note that the "iteration + 100" is a simple way to avoid entering
293 // "if (iteration == 0)" and similar blocks, and also makes it a little
294 // easier to spot the iteration residuals in the DBG file. A more sophisticated
295 // approach can be done later.
296 auto rep = model_.nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver);
297 report += rep;
298 if (rep.converged) {
299 report.converged = true;
300 }
301 return report;
302 }
303
306 {
307 return local_reports_accumulated_;
308 }
309
310private:
312 std::pair<SimulatorReportSingle, ConvergenceReport>
313 solveDomain(const Domain& domain,
315 [[maybe_unused]] const int global_iteration,
316 const bool initial_assembly_required = false)
317 {
318 auto& ebosSimulator = model_.ebosSimulator();
319
321 Dune::Timer solveTimer;
322 solveTimer.start();
323 Dune::Timer detailTimer;
324
325 ebosSimulator.model().newtonMethod().setIterationIndex(0);
326
327 // When called, if assembly has already been performed
328 // with the initial values, we only need to check
329 // for local convergence. Otherwise, we must do a local
330 // assembly.
331 int iter = 0;
333 detailTimer.start();
334 ebosSimulator.model().newtonMethod().setIterationIndex(iter);
335 // TODO: we should have a beginIterationLocal function()
336 // only handling the well model for now
337 ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(),
338 ebosSimulator.timeStepSize(),
339 domain);
340 // Assemble reservoir locally.
341 report += this->assembleReservoirDomain(domain);
342 report.assemble_time += detailTimer.stop();
343 }
344 detailTimer.reset();
345 detailTimer.start();
346 std::vector<double> resnorms;
347 auto convreport = this->getDomainConvergence(domain, timer, 0, resnorms);
348 if (convreport.converged()) {
349 // TODO: set more info, timing etc.
350 report.converged = true;
351 return { report, convreport };
352 }
353
354 // We have already assembled for the first iteration,
355 // but not done the Schur complement for the wells yet.
356 detailTimer.reset();
357 detailTimer.start();
358 model_.wellModel().linearizeDomain(domain,
359 ebosSimulator.model().linearizer().jacobian(),
360 ebosSimulator.model().linearizer().residual());
361 const double tt1 = detailTimer.stop();
362 report.assemble_time += tt1;
363 report.assemble_time_well += tt1;
364
365 // Local Newton loop.
366 const int max_iter = model_.param().max_local_solve_iterations_;
367 const auto& grid = ebosSimulator.vanguard().grid();
368 do {
369 // Solve local linear system.
370 // Note that x has full size, we expect it to be nonzero only for in-domain cells.
371 const int nc = grid.size(0);
372 BVector x(nc);
373 detailTimer.reset();
374 detailTimer.start();
375 this->solveJacobianSystemDomain(domain, x);
376 model_.wellModel().postSolveDomain(x, domain);
377 report.linear_solve_time += detailTimer.stop();
378 report.linear_solve_setup_time += model_.linearSolveSetupTime();
379 report.total_linear_iterations = model_.linearIterationsLastSolve();
380
381 // Update local solution. // TODO: x is still full size, should we optimize it?
382 detailTimer.reset();
383 detailTimer.start();
384 this->updateDomainSolution(domain, x);
385 report.update_time += detailTimer.stop();
386
387 // Assemble well and reservoir.
388 detailTimer.reset();
389 detailTimer.start();
390 ++iter;
391 ebosSimulator.model().newtonMethod().setIterationIndex(iter);
392 // TODO: we should have a beginIterationLocal function()
393 // only handling the well model for now
394 // Assemble reservoir locally.
395 ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(),
396 ebosSimulator.timeStepSize(),
397 domain);
398 report += this->assembleReservoirDomain(domain);
399 report.assemble_time += detailTimer.stop();
400
401 // Check for local convergence.
402 detailTimer.reset();
403 detailTimer.start();
404 convreport = this->getDomainConvergence(domain, timer, iter, resnorms);
405
406 // apply the Schur complement of the well model to the
407 // reservoir linearized equations
408 detailTimer.reset();
409 detailTimer.start();
410 model_.wellModel().linearizeDomain(domain,
411 ebosSimulator.model().linearizer().jacobian(),
412 ebosSimulator.model().linearizer().residual());
413 const double tt2 = detailTimer.stop();
414 report.assemble_time += tt2;
415 report.assemble_time_well += tt2;
416 } while (!convreport.converged() && iter <= max_iter);
417
418 ebosSimulator.problem().endIteration();
419
420 report.converged = convreport.converged();
421 report.total_newton_iterations = iter;
422 report.total_linearizations = iter;
423 report.total_time = solveTimer.stop();
424 // TODO: set more info, timing etc.
425 return { report, convreport };
426 }
427
429 SimulatorReportSingle assembleReservoirDomain(const Domain& domain)
430 {
431 // -------- Mass balance equations --------
432 model_.ebosSimulator().model().linearizer().linearizeDomain(domain);
433 return model_.wellModel().lastReport();
434 }
435
437 void solveJacobianSystemDomain(const Domain& domain, BVector& global_x)
438 {
439 const auto& ebosSimulator = model_.ebosSimulator();
440
441 Dune::Timer perfTimer;
442 perfTimer.start();
443
444 const Mat& main_matrix = ebosSimulator.model().linearizer().jacobian().istlMatrix();
445 if (domain_matrices_[domain.index]) {
446 Details::copySubMatrix(main_matrix, domain.cells, *domain_matrices_[domain.index]);
447 } else {
448 domain_matrices_[domain.index] = std::make_unique<Mat>(Details::extractMatrix(main_matrix, domain.cells));
449 }
450 auto& jac = *domain_matrices_[domain.index];
451 auto res = Details::extractVector(ebosSimulator.model().linearizer().residual(),
452 domain.cells);
453 auto x = res;
454
455 // set initial guess
456 global_x = 0.0;
457 x = 0.0;
458
459 auto& linsolver = domain_linsolvers_[domain.index];
460
461 linsolver.prepare(jac, res);
462 model_.linearSolveSetupTime() = perfTimer.stop();
463 linsolver.setResidual(res);
464 linsolver.solve(x);
465
466 Details::setGlobal(x, domain.cells, global_x);
467 }
468
470 void updateDomainSolution(const Domain& domain, const BVector& dx)
471 {
472 auto& ebosSimulator = model_.ebosSimulator();
473 auto& ebosNewtonMethod = ebosSimulator.model().newtonMethod();
474 SolutionVector& solution = ebosSimulator.model().solution(/*timeIdx=*/0);
475
476 ebosNewtonMethod.update_(/*nextSolution=*/solution,
477 /*curSolution=*/solution,
478 /*update=*/dx,
479 /*resid=*/dx,
480 domain.cells); // the update routines of the black
481 // oil model do not care about the
482 // residual
483
484 // if the solution is updated, the intensive quantities need to be recalculated
485 ebosSimulator.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
486 }
487
489 std::pair<double, double> localDomainConvergenceData(const Domain& domain,
490 std::vector<Scalar>& R_sum,
491 std::vector<Scalar>& maxCoeff,
492 std::vector<Scalar>& B_avg,
493 std::vector<int>& maxCoeffCell)
494 {
495 const auto& ebosSimulator = model_.ebosSimulator();
496
497 double pvSumLocal = 0.0;
498 double numAquiferPvSumLocal = 0.0;
499 const auto& ebosModel = ebosSimulator.model();
500 const auto& ebosProblem = ebosSimulator.problem();
501
502 const auto& ebosResid = ebosSimulator.model().linearizer().residual();
503
504 ElementContext elemCtx(ebosSimulator);
505 const auto& gridView = domain.view;
506 const auto& elemEndIt = gridView.template end</*codim=*/0>();
507 IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
508
509 for (auto elemIt = gridView.template begin</*codim=*/0>();
510 elemIt != elemEndIt;
511 ++elemIt)
512 {
513 if (elemIt->partitionType() != Dune::InteriorEntity) {
514 continue;
515 }
516 const auto& elem = *elemIt;
517 elemCtx.updatePrimaryStencil(elem);
518 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
519
520 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
521 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
522 const auto& fs = intQuants.fluidState();
523
524 const auto pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) *
525 ebosModel.dofTotalVolume(cell_idx);
526 pvSumLocal += pvValue;
527
528 if (isNumericalAquiferCell(elem))
529 {
530 numAquiferPvSumLocal += pvValue;
531 }
532
533 model_.getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue,
534 B_avg, R_sum, maxCoeff, maxCoeffCell);
535 }
536
537 // compute local average in terms of global number of elements
538 const int bSize = B_avg.size();
539 for ( int i = 0; i<bSize; ++i )
540 {
541 B_avg[ i ] /= Scalar(domain.cells.size());
542 }
543
544 return {pvSumLocal, numAquiferPvSumLocal};
545 }
546
547 ConvergenceReport getDomainReservoirConvergence(const double reportTime,
548 const double dt,
549 const int iteration,
550 const Domain& domain,
551 std::vector<Scalar>& B_avg,
552 std::vector<Scalar>& residual_norms)
553 {
554 using Vector = std::vector<Scalar>;
555
556 const int numComp = numEq;
557 Vector R_sum(numComp, 0.0 );
558 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest() );
559 std::vector<int> maxCoeffCell(numComp, -1);
560 const auto [ pvSum, numAquiferPvSum]
561 = this->localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell);
562
563 auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt);
564 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
565
566 const double tol_mb = model_.param().local_tolerance_scaling_mb_ *
567 model_.param().tolerance_mb_;
568 // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
569 // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
570 // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
571 const bool use_relaxed = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
572 iteration >= model_.param().min_strict_cnv_iter_;
573 // Tighter bound for local convergence should increase the
574 // likelyhood of: local convergence => global convergence
575 const double tol_cnv = model_.param().local_tolerance_scaling_cnv_
576 * (use_relaxed ? model_.param().tolerance_cnv_relaxed_
577 : model_.param().tolerance_cnv_);
578
579 // Finish computation
580 std::vector<Scalar> CNV(numComp);
581 std::vector<Scalar> mass_balance_residual(numComp);
582 for (int compIdx = 0; compIdx < numComp; ++compIdx )
583 {
584 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
585 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
586 residual_norms.push_back(CNV[compIdx]);
587 }
588
589 // Create convergence report.
590 ConvergenceReport report{reportTime};
591 using CR = ConvergenceReport;
592 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
593 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
594 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
595 CR::ReservoirFailure::Type::Cnv };
596 double tol[2] = { tol_mb, tol_cnv };
597 for (int ii : {0, 1}) {
598 if (std::isnan(res[ii])) {
599 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
600 if (model_.terminalOutputEnabled()) {
601 OpmLog::debug("NaN residual for " + model_.compNames().name(compIdx) + " equation.");
602 }
603 } else if (res[ii] > model_.param().max_residual_allowed_) {
604 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
605 if (model_.terminalOutputEnabled()) {
606 OpmLog::debug("Too large residual for " + model_.compNames().name(compIdx) + " equation.");
607 }
608 } else if (res[ii] < 0.0) {
609 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
610 if (model_.terminalOutputEnabled()) {
611 OpmLog::debug("Negative residual for " + model_.compNames().name(compIdx) + " equation.");
612 }
613 } else if (res[ii] > tol[ii]) {
614 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
615 }
616 }
617 }
618
619 // Output of residuals.
620 if (model_.terminalOutputEnabled())
621 {
622 // Only rank 0 does print to std::cout
623 if (iteration == 0) {
624 std::string msg = fmt::format("Domain {}, size {}, containing cell {}\n| Iter",
625 domain.index, domain.cells.size(), domain.cells[0]);
626 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
627 msg += " MB(";
628 msg += model_.compNames().name(compIdx)[0];
629 msg += ") ";
630 }
631 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
632 msg += " CNV(";
633 msg += model_.compNames().name(compIdx)[0];
634 msg += ") ";
635 }
636 OpmLog::debug(msg);
637 }
638 std::ostringstream ss;
639 ss << "| ";
640 const std::streamsize oprec = ss.precision(3);
641 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
642 ss << std::setw(4) << iteration;
643 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
644 ss << std::setw(11) << mass_balance_residual[compIdx];
645 }
646 for (int compIdx = 0; compIdx < numComp; ++compIdx) {
647 ss << std::setw(11) << CNV[compIdx];
648 }
649 ss.precision(oprec);
650 ss.flags(oflags);
651 OpmLog::debug(ss.str());
652 }
653
654 return report;
655 }
656
657 ConvergenceReport getDomainConvergence(const Domain& domain,
658 const SimulatorTimerInterface& timer,
659 const int iteration,
660 std::vector<double>& residual_norms)
661 {
662 std::vector<Scalar> B_avg(numEq, 0.0);
663 auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
664 timer.currentStepLength(),
665 iteration,
666 domain,
667 B_avg,
668 residual_norms);
669 report += model_.wellModel().getDomainWellConvergence(domain, B_avg);
670 return report;
671 }
672
674 std::vector<int> getSubdomainOrder()
675 {
676 const auto& ebosSimulator = model_.ebosSimulator();
677 const auto& solution = ebosSimulator.model().solution(0);
678
679 std::vector<int> domain_order(domains_.size());
680 switch (model_.param().local_solve_approach_) {
681 case DomainSolveApproach::GaussSeidel: {
682 switch (model_.param().local_domain_ordering_) {
683 case DomainOrderingMeasure::AveragePressure: {
684 // Use average pressures to order domains.
685 std::vector<std::pair<double, int>> avgpress_per_domain(domains_.size());
686 for (const auto& domain : domains_) {
687 double press_sum = 0.0;
688 for (const int c : domain.cells) {
689 press_sum += solution[c][Indices::pressureSwitchIdx];
690 }
691 const double avgpress = press_sum / domain.cells.size();
692 avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index);
693 }
694 // Lexicographical sort by pressure, then index.
695 std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end());
696 // Reverse since we want high-pressure regions solved first.
697 std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end());
698 for (std::size_t ii = 0; ii < domains_.size(); ++ii) {
699 domain_order[ii] = avgpress_per_domain[ii].second;
700 }
701 break;
702 }
703 case DomainOrderingMeasure::Residual: {
704 // Use maximum residual to order domains.
705 const auto& residual = ebosSimulator.model().linearizer().residual();
706 const int num_vars = residual[0].size();
707 std::vector<std::pair<double, int>> maxres_per_domain(domains_.size());
708 for (const auto& domain : domains_) {
709 double maxres = 0.0;
710 for (const int c : domain.cells) {
711 for (int ii = 0; ii < num_vars; ++ii) {
712 maxres = std::max(maxres, std::fabs(residual[c][ii]));
713 }
714 }
715 maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index);
716 }
717 // Lexicographical sort by pressure, then index.
718 std::sort(maxres_per_domain.begin(), maxres_per_domain.end());
719 // Reverse since we want high-pressure regions solved first.
720 std::reverse(maxres_per_domain.begin(), maxres_per_domain.end());
721 for (std::size_t ii = 0; ii < domains_.size(); ++ii) {
722 domain_order[ii] = maxres_per_domain[ii].second;
723 }
724 }
725 break;
726 }
727 break;
728 }
729
730 case DomainSolveApproach::Jacobi:
731 default:
732 std::iota(domain_order.begin(), domain_order.end(), 0);
733 break;
734 }
735
736 return domain_order;
737 }
738
739 template<class GlobalEqVector>
740 void solveDomainJacobi(GlobalEqVector& solution,
741 GlobalEqVector& locally_solved,
742 SimulatorReportSingle& local_report,
743 const int iteration,
744 const SimulatorTimerInterface& timer,
745 const Domain& domain)
746 {
747 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
748 auto initial_local_solution = Details::extractVector(solution, domain.cells);
749 auto res = solveDomain(domain, timer, iteration);
750 local_report = res.first;
751 if (local_report.converged) {
752 auto local_solution = Details::extractVector(solution, domain.cells);
753 Details::setGlobal(local_solution, domain.cells, locally_solved);
754 Details::setGlobal(initial_local_solution, domain.cells, solution);
755 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
756 } else {
757 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
758 Details::setGlobal(initial_local_solution, domain.cells, solution);
759 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
760 }
761 }
762
763 template<class GlobalEqVector>
764 void solveDomainGaussSeidel(GlobalEqVector& solution,
765 GlobalEqVector& locally_solved,
766 SimulatorReportSingle& local_report,
767 const int iteration,
768 const SimulatorTimerInterface& timer,
769 const Domain& domain)
770 {
771 auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain);
772 auto initial_local_solution = Details::extractVector(solution, domain.cells);
773 auto res = solveDomain(domain, timer, iteration);
774 local_report = res.first;
775 if (!local_report.converged) {
776 // We look at the detailed convergence report to evaluate
777 // if we should accept the unconverged solution.
778 const auto& convrep = res.second;
779 // We do not accept a solution if the wells are unconverged.
780 if (!convrep.wellFailed()) {
781 // Calculare the sums of the mb and cnv failures.
782 double mb_sum = 0.0;
783 double cnv_sum = 0.0;
784 for (const auto& rc : convrep.reservoirConvergence()) {
785 if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
786 mb_sum += rc.value();
787 } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) {
788 cnv_sum += rc.value();
789 }
790 }
791 // If not too high, we overrule the convergence failure.
792 const double acceptable_local_mb_sum = 1e-3;
793 const double acceptable_local_cnv_sum = 1.0;
794 if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) {
795 local_report.converged = true;
796 OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index));
797 }
798 }
799 }
800 if (local_report.converged) {
801 auto local_solution = Details::extractVector(solution, domain.cells);
802 Details::setGlobal(local_solution, domain.cells, locally_solved);
803 } else {
804 model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
805 Details::setGlobal(initial_local_solution, domain.cells, solution);
806 model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain);
807 }
808 }
809
810 double computeCnvErrorPvLocal(const Domain& domain,
811 const std::vector<Scalar>& B_avg, double dt) const
812 {
813 double errorPV{};
814 const auto& ebosSimulator = model_.ebosSimulator();
815 const auto& ebosModel = ebosSimulator.model();
816 const auto& ebosProblem = ebosSimulator.problem();
817 const auto& ebosResid = ebosSimulator.model().linearizer().residual();
818
819 for (const int cell_idx : domain.cells) {
820 const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) *
821 ebosModel.dofTotalVolume(cell_idx);
822 const auto& cellResidual = ebosResid[cell_idx];
823 bool cnvViolated = false;
824
825 for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) {
826 using std::fabs;
827 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
828 cnvViolated = cnvViolated || (fabs(CNV) > model_.param().tolerance_cnv_);
829 }
830
831 if (cnvViolated) {
832 errorPV += pvValue;
833 }
834 }
835 return errorPV;
836 }
837
838 BlackoilModelEbos<TypeTag>& model_;
839 std::vector<Domain> domains_;
840 std::vector<std::unique_ptr<Mat>> domain_matrices_;
841 std::vector<ISTLSolverType> domain_linsolvers_;
842 SimulatorReportSingle local_reports_accumulated_;
843};
844
845} // namespace Opm
846
847#endif // OPM_BLACKOILMODELEBOS_NLDD_HEADER_INCLUDED
Definition AquiferInterface.hpp:35
A NLDD implementation for three-phase black oil.
Definition BlackoilModelEbosNldd.hpp:72
void prepareStep()
Called before starting a time step.
Definition BlackoilModelEbosNldd.hpp:180
SimulatorReportSingle nonlinearIterationNldd(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Do one non-linear NLDD iteration.
Definition BlackoilModelEbosNldd.hpp:188
BlackoilModelEbosNldd(BlackoilModelEbos< TypeTag > &model)
The constructor sets up the subdomains.
Definition BlackoilModelEbosNldd.hpp:93
const SimulatorReportSingle & localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModelEbosNldd.hpp:305
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 BlackoilPhases.hpp:27
std::pair< std::vector< int >, int > partitionCells(const Grid &grid, const std::vector< Well > &wells, const std::string &method, const int num_local_domains, const double partition_imbalance)
Partitions the grid using the specified method.
Definition partitionCells.cpp:77
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:237
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34