Skip to content

Commit

Permalink
Merge pull request #1829 from GitPaean/different_control_different_to…
Browse files Browse the repository at this point in the history 
…lerance

letting getResidualMeasureValue() use the same criterion as getWellConvergence()
  • Loading branch information
@atgeirr
atgeirr authored May 16, 2019
2 parents 41ef80d + 080e840 commit 4e5aab5
Show file tree
Hide file tree
Showing 6 changed files with 188 additions and 97 deletions.
8 changes: 7 additions & 1 deletion opm/simulators/wells/MultisegmentWell.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,7 +379,13 @@ namespace Opm
void detectOscillations(const std::vector<double>& measure_history,
const int it, bool& oscillate, bool& stagnate) const;

double getResidualMeasureValue(const std::vector<double>& residuals) const;
double getResidualMeasureValue(const std::vector<double>& residuals,
DeferredLogger& deferred_logger) const;

double getControlTolerance(DeferredLogger& deferred_logger) const;

void checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const;

};

Expand Down
117 changes: 82 additions & 35 deletions opm/simulators/wells/MultisegmentWell_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -429,28 +429,9 @@ namespace Opm
}
}

using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
switch(well_controls_get_current_type(well_controls_)) {
case THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
break;
case BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
ctrltype = CR::WellFailure::Type::ControlRate;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
}
assert(ctrltype != CR::WellFailure::Type::Invalid);

std::vector<double> maximum_residual(numWellEq, 0.0);

ConvergenceReport report;
const int dummy_component = -1;
// TODO: the following is a little complicated, maybe can be simplified in some way?
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
for (int seg = 0; seg < numberOfSegments(); ++seg) {
Expand All @@ -460,18 +441,8 @@ namespace Opm
maximum_residual[eq_idx] = flux_residual;
}
} else { // pressure or control equation
if (seg == 0) {
// Control equation
const double control_residual = abs_residual[seg][eq_idx];
if (std::isnan(control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
} else if (control_residual > param_.max_residual_allowed_) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
// TODO: we should distinguish the flux residual or pressure residual here
} else if (control_residual > param_.tolerance_wells_) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
}
} else {
// for the top segment (seg == 0), it is control equation, will be checked later separately
if (seg > 0) {
// Pressure equation
const double pressure_residual = abs_residual[seg][eq_idx];
if (pressure_residual > maximum_residual[eq_idx]) {
Expand All @@ -482,6 +453,7 @@ namespace Opm
}
}

using CR = ConvergenceReport;
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
if (eq_idx < num_components_) { // phase or component mass equations
const double flux_residual = maximum_residual[eq_idx];
Expand All @@ -495,6 +467,7 @@ namespace Opm
}
} else { // pressure equation
const double pressure_residual = maximum_residual[eq_idx];
const int dummy_component = -1;
if (std::isnan(pressure_residual)) {
report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()});
} else if (std::isinf(pressure_residual)) {
Expand All @@ -505,6 +478,8 @@ namespace Opm
}
}

checkConvergenceControlEq(report, deferred_logger);

return report;
}

Expand Down Expand Up @@ -1825,7 +1800,7 @@ namespace Opm
}

residual_history.push_back(getWellResiduals(B_avg));
measure_history.push_back(getResidualMeasureValue(residual_history[it]));
measure_history.push_back(getResidualMeasureValue(residual_history[it], deferred_logger) );

bool is_oscillate = false;
bool is_stagnate = false;
Expand Down Expand Up @@ -2176,7 +2151,6 @@ namespace Opm
assert(int(B_avg.size() ) == num_components_);
std::vector<Scalar> residuals(numWellEq + 1, 0.0);

// TODO: maybe we should distinguish the bhp control or rate control equations here
for (int seg = 0; seg < numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
double residual = 0.;
Expand Down Expand Up @@ -2248,7 +2222,8 @@ namespace Opm
template<typename TypeTag>
double
MultisegmentWell<TypeTag>::
getResidualMeasureValue(const std::vector<double>& residuals) const
getResidualMeasureValue(const std::vector<double>& residuals,
DeferredLogger& deferred_logger) const
{
assert(int(residuals.size()) == numWellEq + 1);

Expand All @@ -2268,11 +2243,83 @@ namespace Opm
++count;
}

const double control_tolerance = getControlTolerance(deferred_logger);
if (residuals[SPres + 1] > control_tolerance) {
sum += residuals[SPres + 1] / control_tolerance;
++count;
}

// if (count == 0), it should be converged.
assert(count != 0);

// return sum / double(count);
return sum;
}





template<typename TypeTag>
double
MultisegmentWell<TypeTag>::
getControlTolerance(DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
switch(well_controls_get_current_type(well_controls_) ) {
case BHP:
control_tolerance = param_.tolerance_wells_;
break;
case THP:
control_tolerance = param_.tolerance_pressure_ms_wells_;
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
control_tolerance = param_.tolerance_wells_;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
}

return control_tolerance;
}




template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const
{
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
switch(well_controls_get_current_type(well_controls_)) {
case THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
break;
case BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
ctrltype = CR::WellFailure::Type::ControlRate;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
}
assert(ctrltype != CR::WellFailure::Type::Invalid);

const double control_residual = std::abs(resWell_[0][SPres]);
const double control_tolerance = getControlTolerance(deferred_logger);

const int dummy_component = -1;
if (std::isnan(control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
} else if (control_residual > param_.max_residual_allowed_) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
} else if (control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
}
}
}
4 changes: 4 additions & 0 deletions opm/simulators/wells/StandardWell.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -439,6 +439,10 @@ namespace Opm
Opm::DeferredLogger& deferred_logger) override;

virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;

void checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const;

};

}
Expand Down
3 changes: 3 additions & 0 deletions opm/simulators/wells/StandardWellV.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -465,6 +465,9 @@ namespace Opm
Opm::DeferredLogger& deferred_logger);

virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;

void checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
};

}
Expand Down
77 changes: 47 additions & 30 deletions opm/simulators/wells/StandardWellV_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2029,40 +2029,12 @@ namespace Opm
}
}

// processing the residual of the well control equation
const double well_control_residual = res[numWellEq_ - 1];
// TODO: we should have better way to specify the control equation tolerance
double control_tolerance = 0.;
switch(well_controls_get_current_type(well_controls_)) {
case THP:
type = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case BHP: // pressure type of control
type = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
type = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; // smaller tolerance for rate control
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " + name(), deferred_logger);
}

const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({type, CR::Severity::NotANumber, dummy_component, name()});
} else if (well_control_residual > maxResidualAllowed * 10.) {
report.setWellFailed({type, CR::Severity::TooLarge, dummy_component, name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
}
checkConvergenceControlEq(report, deferred_logger);

if (this->has_polymermw && well_type_ == INJECTOR) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
Expand Down Expand Up @@ -3112,4 +3084,49 @@ namespace Opm
duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
}
}





template<typename TypeTag>
void
StandardWellV<TypeTag>::
checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
switch(well_controls_get_current_type(well_controls_)) {
case THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case BHP: // pressure type of control
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; // smaller tolerance for rate control
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
}
assert(ctrltype != CR::WellFailure::Type::Invalid);

const double well_control_residual = std::abs(resWell_[0][Bhp]);
const int dummy_component = -1;
const double max_residual_allowed = param_.max_residual_allowed_;
if (std::isnan(well_control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
} else if (well_control_residual > max_residual_allowed * 10.) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
}

}
}
Loading

0 comments on commit 4e5aab5

Please sign in to comment.