OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4739 from akva2/blackoilmodelebos_nldd_solver_approach

Browse Source 
BlackoilModelEbos: improve nldd solver approach
This commit is contained in:
Bård Skaflestad 2023-07-03 17:40:45 +02:00 committed by GitHub
commit 4498d14a78
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 106 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@ -541,60 +541,16 @@ namespace Opm {
for (const int domain_index : domain_order) { for (const int domain_index : domain_order) {
const auto& domain = domains_[domain_index]; const auto& domain = domains_[domain_index];
SimulatorReportSingle local_report; SimulatorReportSingle local_report;
if (param_.local_solve_approach_ == "jacobi") { switch (param_.local_solve_approach_) {
auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); case DomainSolveApproach::Jacobi:
auto initial_local_solution = Details::extractVector(solution, domain.cells); solveDomainJacobi(solution, locally_solved, local_report,
auto res = solveDomain(domain, timer, iteration); iteration, timer, domain);
local_report = res.first; break;
if (local_report.converged) { default:
auto local_solution = Details::extractVector(solution, domain.cells); case DomainSolveApproach::GaussSeidel:
Details::setGlobal(local_solution, domain.cells, locally_solved); solveDomainGaussSeidel(solution, locally_solved, local_report,
Details::setGlobal(initial_local_solution, domain.cells, solution); iteration, timer, domain);
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); break;
} else {
wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
Details::setGlobal(initial_local_solution, domain.cells, solution);
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view);
}
} 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);
}
} }
// This should have updated the global matrix to be // This should have updated the global matrix to be
// dR_i/du_j evaluated at new local solutions for // dR_i/du_j evaluated at new local solutions for
@ -624,7 +580,7 @@ namespace Opm {
local_reports_accumulated_ += rep; local_reports_accumulated_ += rep;
} }
if (param_.local_solve_approach_ == "jacobi") { if (param_.local_solve_approach_ == DomainSolveApproach::Jacobi) {
solution = locally_solved; solution = locally_solved;
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
} }
@ -1794,7 +1750,8 @@ namespace Opm {
const auto& solution = ebosSimulator().model().solution(0); const auto& solution = ebosSimulator().model().solution(0);
std::vector<int> domain_order(domains_.size()); std::vector<int> domain_order(domains_.size());
if (param_.local_solve_approach_ == "gauss-seidel") { switch (param_.local_solve_approach_) {
case DomainSolveApproach::GaussSeidel: {
switch (param_.local_domain_ordering_) { switch (param_.local_domain_ordering_) {
case DomainOrderingMeasure::AveragePressure: { case DomainOrderingMeasure::AveragePressure: {
// Use average pressures to order domains. // Use average pressures to order domains.
@ -1838,14 +1795,91 @@ namespace Opm {
domain_order[ii] = maxres_per_domain[ii].second; domain_order[ii] = maxres_per_domain[ii].second;
} }
} }
break;
} }
} else { break;
}
case DomainSolveApproach::Jacobi:
default:
std::iota(domain_order.begin(), domain_order.end(), 0); std::iota(domain_order.begin(), domain_order.end(), 0);
break;
} }
return domain_order; return domain_order;
} }
template<class GlobalEqVector>
void solveDomainJacobi(GlobalEqVector& solution,
GlobalEqVector& locally_solved,
SimulatorReportSingle& local_report,
const int iteration,
const SimulatorTimerInterface& timer,
const Domain& domain)
{
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) {
auto local_solution = Details::extractVector(solution, domain.cells);
Details::setGlobal(local_solution, domain.cells, locally_solved);
Details::setGlobal(initial_local_solution, domain.cells, solution);
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view);
} else {
wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars);
Details::setGlobal(initial_local_solution, domain.cells, solution);
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view);
}
}
template<class GlobalEqVector>
void solveDomainGaussSeidel(GlobalEqVector& solution,
GlobalEqVector& locally_solved,
SimulatorReportSingle& local_report,
const int iteration,
const SimulatorTimerInterface& timer,
const Domain& domain)
{
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: public:
std::vector<bool> wasSwitched_; std::vector<bool> wasSwitched_;
}; };

12
opm/simulators/flow/BlackoilModelParametersEbos.hpp
View File

@ -539,7 +539,7 @@ namespace Opm
/// Nonlinear solver type: newton or nldd. /// Nonlinear solver type: newton or nldd.
std::string nonlinear_solver_; std::string nonlinear_solver_;
/// 'jacobi' and 'gauss-seidel' supported. /// 'jacobi' and 'gauss-seidel' supported.
std::string local_solve_approach_; DomainSolveApproach local_solve_approach_{DomainSolveApproach::Jacobi};
int max_local_solve_iterations_; int max_local_solve_iterations_;
@ -587,7 +587,15 @@ namespace Opm
max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches); max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells); use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver); nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver);
local_solve_approach_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach); std::string approach = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach);
if (approach == "jacobi") {
local_solve_approach_ = DomainSolveApproach::Jacobi;
} else if (approach == "gauss-seidel") {
local_solve_approach_ = DomainSolveApproach::GaussSeidel;
} else {
throw std::runtime_error("Invalid domain solver approach '" + approach + "' specified.");
}
max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations); max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations);
local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb); local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb);
local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv); local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv);

5
opm/simulators/flow/SubDomain.hpp
View File

@ -27,6 +27,11 @@
namespace Opm namespace Opm
{ {
//! \brief Solver approach for NLDD.
enum class DomainSolveApproach {
Jacobi,
GaussSeidel
};
//! \brief Measure to use for domain ordering. //! \brief Measure to use for domain ordering.
enum class DomainOrderingMeasure { enum class DomainOrderingMeasure {