Skip to content

Commit

Permalink
Report CNV Violation Pore-Volume Fraction to INFOITER
Browse files Browse the repository at this point in the history 
This commit includes the fraction of pore-volume whose CNV targets
are violated as a new per-iteration quantity in the INFOITER file
(--output-extra-convergence-info=iteration), with the column header
"CnvErrPvFrac".  We collect the values which are already calculated
in

    BlackoilModel<>::getReservoirConvergence()

and store these as a pair of numerator and denominator in the
ConvergenceReport class.  Note that we need both the numerator and
the denominator in order to aggregate contributions from multiple
ranks.

While here, also make a few more objects 'const' and calculate
column widths directly instead of the maximum number of characters
in writeConvergenceHeader().
  • Loading branch information
@bska
bska committed Apr 29, 2024
1 parent 5bdb1d4 commit 4d646f8
Show file tree
Hide file tree
Showing 7 changed files with 190 additions and 61 deletions.
22 changes: 13 additions & 9 deletions opm/simulators/flow/BlackoilModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -897,17 +897,18 @@ namespace Opm {
Vector R_sum(numComp, 0.0 );
Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
std::vector<int> maxCoeffCell(numComp, -1);
const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
const auto [ pvSumLocal, numAquiferPvSumLocal] =
this->localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);

// compute global sum and max of quantities
const auto [ pvSum, numAquiferPvSum ] =
convergenceReduction(grid_.comm(), pvSumLocal,
numAquiferPvSumLocal,
R_sum, maxCoeff, B_avg);

auto cnvErrorPvFraction = computeCnvErrorPv(B_avg, dt);
cnvErrorPvFraction /= (pvSum - numAquiferPvSum);

const auto cnvErrorPvFraction =
computeCnvErrorPv(B_avg, dt)
/ (pvSum - numAquiferPvSum);

// For each iteration, we need to determine whether to use the relaxed tolerances.
// To disable the usage of relaxed tolerances, you can set the relaxed tolerances as the strict tolerances.
Expand All @@ -927,7 +928,7 @@ namespace Opm {
const bool relax_pv_fraction_cnv = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_;
const bool use_relaxed_cnv = relax_final_iteration_cnv || relax_pv_fraction_cnv || relax_iter_cnv;

if (relax_final_iteration_mb || relax_final_iteration_cnv) {
if (relax_final_iteration_mb || relax_final_iteration_cnv) {
if ( terminal_output_ ) {
std::string message = "Number of newton iterations reached its maximum try to continue with relaxed tolerances:";
if (relax_final_iteration_mb)
Expand All @@ -944,7 +945,7 @@ namespace Opm {
// Finish computation
std::vector<Scalar> CNV(numComp);
std::vector<Scalar> mass_balance_residual(numComp);
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
for (int compIdx = 0; compIdx < numComp; ++compIdx)
{
CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
Expand All @@ -953,12 +954,15 @@ namespace Opm {

// Create convergence report.
ConvergenceReport report{reportTime};
report.setPoreVolCnvViolationFraction(cnvErrorPvFraction, pvSum - numAquiferPvSum);

using CR = ConvergenceReport;
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
const double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
const CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
CR::ReservoirFailure::Type::Cnv };
double tol[2] = { tol_mb, tol_cnv };
const double tol[2] = { tol_mb, tol_cnv };

for (int ii : {0, 1}) {
if (std::isnan(res[ii])) {
report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
Expand Down
1 change: 1 addition & 0 deletions opm/simulators/flow/BlackoilModelNldd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -656,6 +656,7 @@ class BlackoilModelNldd {
} else if (res[ii] > tol[ii]) {
report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
}

report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
}
}
Expand Down
83 changes: 63 additions & 20 deletions opm/simulators/flow/ExtraConvergenceOutputThread.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>

#include <algorithm>
#include <array>
#include <cassert>
#include <condition_variable>
#include <cstddef>
Expand All @@ -40,12 +41,44 @@
#include <stdexcept>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

auto fixedHeaders() noexcept
{
using namespace std::literals::string_literals;

return std::array {
"ReportStep"s,
"TimeStep"s,
"Time"s,
"CnvErrPvFrac"s,
"Iteration"s,
};
}

template <typename HeaderSequence>
auto maxHeaderSize(const HeaderSequence& headers)
{
using sz_t = std::remove_cv_t<std::remove_reference_t<
decltype(std::declval<HeaderSequence>().front().size())>>;

if (headers.empty()) {
return sz_t{0};
}

return std::accumulate(headers.begin(), headers.end(), sz_t{1},
[](const auto m, const auto& header)
{
return std::max(m, header.size() + 1);
});
}

std::string
formatMetricColumn(const Opm::ConvergenceOutputThread::ComponentToPhaseName& getPhaseName,
const Opm::ConvergenceReport::ReservoirConvergenceMetric& metric)
Expand All @@ -65,45 +98,49 @@ namespace {
[&getPhaseName](const std::string::size_type maxChar,
const Opm::ConvergenceReport::ReservoirConvergenceMetric& metric)
{
return std::max(maxChar, formatMetricColumn(getPhaseName, metric).size());
return std::max(maxChar, formatMetricColumn(getPhaseName, metric).size() + 1);
});
}

std::string::size_type
std::pair<std::string::size_type, std::string::size_type>
writeConvergenceHeader(std::ostream& os,
const Opm::ConvergenceOutputThread::ComponentToPhaseName& getPhaseName,
const Opm::ConvergenceReportQueue::OutputRequest& firstRequest)
{
const auto minColSize = std::string::size_type{11};
const auto initial_headers = fixedHeaders();
const auto minColSize = maxHeaderSize(initial_headers);

os << std::right << std::setw(minColSize) << "ReportStep" << ' '
<< std::right << std::setw(minColSize) << "TimeStep" << ' '
<< std::right << std::setw(minColSize) << "Time" << ' '
<< std::right << std::setw(minColSize) << "Iteration";
{
auto leadingSpace = false;

const auto& metrics = firstRequest.reports.front().reservoirConvergence();
const auto maxChar = maxColHeaderSize(minColSize, getPhaseName, metrics);
for (const auto& columnHeader : initial_headers) {
if (leadingSpace) { os << ' '; }
os << std::right << std::setw(minColSize) << columnHeader;
leadingSpace = true;
}
}

const auto& metrics = firstRequest.reports.front().reservoirConvergence();
const auto headerSize = maxColHeaderSize(minColSize, getPhaseName, metrics);

for (const auto& metric : metrics) {
os << std::right << std::setw(maxChar + 1)
os << std::right << std::setw(headerSize)
<< formatMetricColumn(getPhaseName, metric);
}

// Note: Newline character intentionally placed in separate output
// request to not influence right-justification of column header.
os << std::right << std::setw(maxChar + 1) << "WellStatus" << '\n';
os << std::right << std::setw(headerSize) << "WellStatus" << '\n';

return maxChar;
return { minColSize, headerSize };
}


void writeConvergenceRequest(std::ostream& os,
const Opm::ConvergenceOutputThread::ConvertToTimeUnits& convertTime,
std::string::size_type colSize,
const std::string::size_type firstColSize,
const std::string::size_type colSize,
const Opm::ConvergenceReportQueue::OutputRequest& request)
{
const auto firstColSize = std::string::size_type{11};

os.setf(std::ios_base::scientific);

auto iter = 0;
Expand All @@ -112,19 +149,23 @@ namespace {
<< std::setw(firstColSize) << request.currentStep << ' '
<< std::setprecision(4) << std::setw(firstColSize)
<< convertTime(report.reportTime()) << ' '
<< std::setprecision(4) << std::setw(firstColSize)
<< report.cnvViolatedPvFraction() << ' '
<< std::setw(firstColSize) << iter;

for (const auto& metric : report.reservoirConvergence()) {
os << std::setprecision(4) << std::setw(colSize + 1) << metric.value();
os << std::setprecision(4) << std::setw(colSize) << metric.value();
}

os << std::right << std::setw(colSize + 1)
os << std::right << std::setw(colSize)
<< (report.wellFailed() ? "FAIL" : "CONV");

if (report.wellFailed()) {
for (const auto& wf : report.wellFailures()) {
os << " " << to_string(wf);
os << ' ' << to_string(wf);
}
}

os << '\n';

++iter;
Expand Down Expand Up @@ -160,6 +201,7 @@ class Opm::ConvergenceOutputThread::Impl
ComponentToPhaseName getPhaseName_{};
ConvertToTimeUnits convertTime_{};
std::optional<std::ofstream> infoIter_{};
std::string::size_type firstColSize_{0};
std::string::size_type colSize_{0};
bool haveOutputIterHeader_{false};
bool finalRequestWritten_{false};
Expand Down Expand Up @@ -203,7 +245,7 @@ writeIterInfo(const std::vector<ConvergenceReportQueue::OutputRequest>& requests
}

if (! this->haveOutputIterHeader_) {
this->colSize_ =
std::tie(this->firstColSize_, this->colSize_) =
writeConvergenceHeader(this->infoIter_.value(),
this->getPhaseName_,
requests.front());
Expand All @@ -213,6 +255,7 @@ writeIterInfo(const std::vector<ConvergenceReportQueue::OutputRequest>& requests
for (const auto& request : requests) {
writeConvergenceRequest(this->infoIter_.value(),
this->convertTime_,
this->firstColSize_,
this->colSize_,
request);

Expand Down
22 changes: 16 additions & 6 deletions opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -363,13 +363,16 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
+ schedule().seconds(timer.currentStepNum()),
timer.currentStepLength());
simulator_.setEpisodeIndex(timer.currentStepNum());

if (serializer_.shouldLoad()) {
wellModel_().prepareDeserialize(serializer_.loadStep() - 1);
serializer_.loadState();
simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
solver_->model().beginReportStep();
bool enableTUNING = Parameters::get<TypeTag, Properties::EnableTuning>();

this->solver_->model().beginReportStep();

const bool enableTUNING = Parameters::get<TypeTag, Properties::EnableTuning>();

// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are too large for the solver to converge
Expand Down Expand Up @@ -442,10 +445,17 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
report_.success.solver_time += solverTimer_->secsSinceStart();

if (this->grid().comm().rank() == 0) {
// Grab the step convergence reports that are new since last we were here.
const auto& reps = solver_->model().stepReports();
this->writeConvergenceOutput(std::vector<StepReport>{reps.begin() + already_reported_steps_, reps.end()});
already_reported_steps_ = reps.size();
// Grab the step convergence reports that are new since last we
// were here.
const auto& reps = this->solver_->model().stepReports();

auto reports = std::vector<StepReport> {
reps.begin() + this->already_reported_steps_, reps.end()
};

this->writeConvergenceOutput(std::move(reports));

this->already_reported_steps_ = reps.size();
}

// Increment timer, remember well state.
Expand Down

0 comments on commit 4d646f8

Please sign in to comment.