diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index b0f832eb3..47a30e3b7 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -545,44 +545,8 @@ namespace Opm { solveDomainJacobi(solution, locally_solved, local_report, iteration, timer, domain); } else { - assert(param_.local_solve_approach_ == "gauss-seidel"); - auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); - auto initial_local_solution = Details::extractVector(solution, domain.cells); - auto res = solveDomain(domain, timer, iteration); - local_report = res.first; - if (!local_report.converged) { - // We look at the detailed convergence report to evaluate - // if we should accept the unconverged solution. - const auto& convrep = res.second; - // We do not accept a solution if the wells are unconverged. - if (!convrep.wellFailed()) { - // Calculare the sums of the mb and cnv failures. - double mb_sum = 0.0; - double cnv_sum = 0.0; - for (const auto& rc : convrep.reservoirConvergence()) { - if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { - mb_sum += rc.value(); - } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { - cnv_sum += rc.value(); - } - } - // If not too high, we overrule the convergence failure. - const double acceptable_local_mb_sum = 1e-3; - const double acceptable_local_cnv_sum = 1.0; - if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) { - local_report.converged = true; - OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index)); - } - } - } - if (local_report.converged) { - auto local_solution = Details::extractVector(solution, domain.cells); - Details::setGlobal(local_solution, domain.cells, locally_solved); - } else { - wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); - Details::setGlobal(initial_local_solution, domain.cells, solution); - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); - } + solveDomainGaussSeidel(solution, locally_solved, local_report, + iteration, timer, domain); } // This should have updated the global matrix to be // dR_i/du_j evaluated at new local solutions for @@ -1858,6 +1822,54 @@ namespace Opm { } } + template + void solveDomainGaussSeidel(GlobalEqVector& solution, + GlobalEqVector& locally_solved, + SimulatorReportSingle& local_report, + const int iteration, + const SimulatorTimerInterface& timer, + const Domain& domain) + { + assert(param_.local_solve_approach_ == "gauss-seidel"); + auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveDomain(domain, timer, iteration); + local_report = res.first; + if (!local_report.converged) { + // We look at the detailed convergence report to evaluate + // if we should accept the unconverged solution. + const auto& convrep = res.second; + // We do not accept a solution if the wells are unconverged. + if (!convrep.wellFailed()) { + // Calculare the sums of the mb and cnv failures. + double mb_sum = 0.0; + double cnv_sum = 0.0; + for (const auto& rc : convrep.reservoirConvergence()) { + if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { + mb_sum += rc.value(); + } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { + cnv_sum += rc.value(); + } + } + // If not too high, we overrule the convergence failure. + const double acceptable_local_mb_sum = 1e-3; + const double acceptable_local_cnv_sum = 1.0; + if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) { + local_report.converged = true; + OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index)); + } + } + } + if (local_report.converged) { + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + } else { + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + } + public: std::vector wasSwitched_; };