BlackoilModelEbos: put Gauss-Seidel NLDD solve in separate method

increases readability of nonlinearIterationNldd
This commit is contained in:
Arne Morten Kvarving 2023-07-03 12:40:15 +02:00
parent faa4bb25d7
commit 8eb5e173c3

View File

@ -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<class GlobalEqVector>
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<bool> wasSwitched_;
};