mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-01 03:37:15 -06:00
582 lines
15 KiB
C++
582 lines
15 KiB
C++
/*
|
|
Copyright 2017 SINTEF ICT, Applied Mathematics.
|
|
Copyright 2017 Statoil ASA.
|
|
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include <examples/exampleSetup.hpp>
|
|
|
|
#include <opm/utility/ECLCaseUtilities.hpp>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cassert>
|
|
#include <cmath>
|
|
#include <cstdlib>
|
|
#include <exception>
|
|
#include <functional>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <map>
|
|
#include <memory>
|
|
#include <numeric>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <type_traits>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include <boost/algorithm/string/case_conv.hpp>
|
|
#include <boost/filesystem.hpp>
|
|
#include <boost/filesystem/fstream.hpp>
|
|
#include <regex>
|
|
|
|
#include <ert/ecl/ecl_file.h>
|
|
#include <ert/ecl/ecl_file_kw.h>
|
|
#include <ert/ecl/ecl_file_view.h>
|
|
#include <ert/ecl/ecl_kw.h>
|
|
#include <ert/ecl/ecl_kw_magic.h>
|
|
#include <ert/util/ert_unique_ptr.hpp>
|
|
|
|
namespace StringUtils {
|
|
namespace {
|
|
std::string trim(const std::string& s)
|
|
{
|
|
const auto anchor_ws =
|
|
std::regex(R"~~(^\s+([^\s]+)\s+$)~~");
|
|
|
|
auto m = std::smatch{};
|
|
|
|
if (std::regex_match(s, m, anchor_ws)) {
|
|
return m[1];
|
|
}
|
|
|
|
return s;
|
|
}
|
|
|
|
std::vector<std::string> split(const std::string& s)
|
|
{
|
|
if (s.empty()) {
|
|
// Single element vector whose only element is the empty
|
|
// string.
|
|
return { "" };
|
|
}
|
|
|
|
const auto sep = std::regex(R"~~([\s,;.|]+)~~");
|
|
|
|
using TI = std::sregex_token_iterator;
|
|
|
|
// vector<string>(begin, end)
|
|
//
|
|
// Range is every substring (i.e., token) in input string 's'
|
|
// that does NOT match 'sep'.
|
|
return{ TI(s.begin(), s.end(), sep, -1), TI{} };
|
|
}
|
|
} // namespace Anonymous
|
|
|
|
template <typename T>
|
|
struct StringTo;
|
|
|
|
template <>
|
|
struct StringTo<int>
|
|
{
|
|
static int value(const std::string& s);
|
|
};
|
|
|
|
template <>
|
|
struct StringTo<double>
|
|
{
|
|
static double value(const std::string& s);
|
|
};
|
|
|
|
template <>
|
|
struct StringTo<std::string>
|
|
{
|
|
static std::string value(const std::string& s);
|
|
};
|
|
|
|
int StringTo<int>::value(const std::string& s)
|
|
{
|
|
return std::stoi(s);
|
|
}
|
|
|
|
double StringTo<double>::value(const std::string& s)
|
|
{
|
|
return std::stod(s);
|
|
}
|
|
|
|
std::string StringTo<std::string>::value(const std::string& s)
|
|
{
|
|
return trim(s);
|
|
}
|
|
|
|
namespace VectorValue {
|
|
template <typename T>
|
|
std::vector<T> get(const std::string& s, std::true_type)
|
|
{
|
|
return split(s);
|
|
}
|
|
|
|
template <typename T>
|
|
std::vector<T> get(const std::string& s, std::false_type)
|
|
{
|
|
const auto tokens = split(s);
|
|
|
|
auto ret = std::vector<T>{};
|
|
ret.reserve(tokens.size());
|
|
for (const auto& token : tokens) {
|
|
ret.push_back(StringTo<T>::value(token));
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
template <typename T>
|
|
std::vector<T> get(const std::string& s)
|
|
{
|
|
return get<T>(s, typename std::is_same<T, std::string>::type());
|
|
}
|
|
} // namespace VectorValue
|
|
} // namespace StringUtils
|
|
|
|
namespace {
|
|
struct Reference
|
|
{
|
|
std::vector<double> data;
|
|
};
|
|
|
|
struct Calculated
|
|
{
|
|
std::vector<double> data;
|
|
};
|
|
|
|
class VectorUnits
|
|
{
|
|
private:
|
|
using USys = ::Opm::ECLUnits::UnitSystem;
|
|
|
|
public:
|
|
using UnitConvention = ::Opm::ECLGraph::UnitConvention;
|
|
|
|
VectorUnits()
|
|
: units_({ { "pressure", &USys::pressure } })
|
|
{
|
|
}
|
|
|
|
UnitConvention getUnit(const std::string& vector) const
|
|
{
|
|
auto p = units_.find(vector);
|
|
|
|
if (p == units_.end()) {
|
|
std::ostringstream os;
|
|
|
|
os << "Unsupported Vector Quantity '" << vector << '\'';
|
|
|
|
throw std::domain_error(os.str());
|
|
}
|
|
|
|
return p->second;
|
|
}
|
|
|
|
private:
|
|
std::map<std::string, UnitConvention> units_;
|
|
};
|
|
|
|
class VectorDifference
|
|
{
|
|
public:
|
|
using Vector = std::vector<double>;
|
|
using size_type = Vector::size_type;
|
|
|
|
VectorDifference(const Vector& x, const Vector& y)
|
|
: x_(x), y_(y)
|
|
{
|
|
if (x_.size() != y_.size()) {
|
|
std::ostringstream os;
|
|
|
|
os << "Incompatible Array Sizes: Expected 2x"
|
|
<< x_.size() << ", but got ("
|
|
<< x_.size() << ", " << y_.size() << ')';
|
|
|
|
throw std::domain_error(os.str());
|
|
}
|
|
}
|
|
|
|
size_type size() const
|
|
{
|
|
return x_.size();
|
|
}
|
|
|
|
bool empty() const
|
|
{
|
|
return this->size() == 0;
|
|
}
|
|
|
|
double operator[](const size_type i) const
|
|
{
|
|
return x_[i] - y_[i];
|
|
}
|
|
|
|
private:
|
|
const Vector& x_;
|
|
const Vector& y_;
|
|
};
|
|
|
|
template <class Vector1, class Vector2>
|
|
class VectorRatio
|
|
{
|
|
public:
|
|
using size_type = typename std::decay<
|
|
decltype(std::declval<Vector1>()[0])
|
|
>::type;
|
|
|
|
VectorRatio(const Vector1& x, const Vector2& y)
|
|
: x_(x), y_(y)
|
|
{
|
|
if (x_.size() != y.size()) {
|
|
std::ostringstream os;
|
|
|
|
os << "Incompatible Array Sizes: Expected 2x"
|
|
<< x_.size() << ", but got ("
|
|
<< x_.size() << ", " << y_.size() << ')';
|
|
|
|
throw std::domain_error(os.str());
|
|
}
|
|
}
|
|
|
|
size_type size() const
|
|
{
|
|
return x_.size();
|
|
}
|
|
|
|
bool empty() const
|
|
{
|
|
return x_.empty();
|
|
}
|
|
|
|
double operator[](const size_type i) const
|
|
{
|
|
return x_[i] / y_[i];
|
|
}
|
|
|
|
private:
|
|
const Vector1& x_;
|
|
const Vector2& y_;
|
|
};
|
|
|
|
struct ErrorMeasurement
|
|
{
|
|
double volume;
|
|
double inf;
|
|
};
|
|
|
|
struct ErrorTolerance
|
|
{
|
|
double absolute;
|
|
double relative;
|
|
};
|
|
|
|
struct AggregateErrors
|
|
{
|
|
std::vector<ErrorMeasurement> absolute;
|
|
std::vector<ErrorMeasurement> relative;
|
|
};
|
|
|
|
struct ReferenceSolution
|
|
{
|
|
std::vector<double> raw;
|
|
std::vector<double> SI;
|
|
};
|
|
|
|
template <class FieldVariable>
|
|
double volumeMetric(const FieldVariable& x)
|
|
{
|
|
auto result = 0.0;
|
|
|
|
for (decltype(x.size())
|
|
i = 0, n = x.size(); i < n; ++i)
|
|
{
|
|
const auto m = std::abs(x[i]);
|
|
result += m * m;
|
|
}
|
|
|
|
return std::sqrt(result / x.size());
|
|
}
|
|
|
|
template <class FieldVariable>
|
|
double pointMetric(const FieldVariable& x)
|
|
{
|
|
static_assert(std::is_same<typename std::decay<decltype(std::abs(x[0]))>::type, double>::value,
|
|
"Field Variable Value Type Must be 'double'");
|
|
|
|
if (x.empty()) {
|
|
return 0;
|
|
}
|
|
|
|
auto max = 0*x[0] - 1;
|
|
|
|
for (decltype(x.size())
|
|
i = 0, n = x.size(); i < n; ++i)
|
|
{
|
|
const auto t = std::abs(x[i]);
|
|
|
|
if (t > max) {
|
|
max = t;
|
|
}
|
|
}
|
|
|
|
return max;
|
|
}
|
|
|
|
ErrorTolerance
|
|
testTolerances(const ::Opm::ParameterGroup& param)
|
|
{
|
|
const auto atol = param.getDefault("atol", 1.0e-8);
|
|
const auto rtol = param.getDefault("rtol", 5.0e-12);
|
|
|
|
return ErrorTolerance{ atol, rtol };
|
|
}
|
|
|
|
std::vector<std::string>
|
|
testQuantities(const ::Opm::ParameterGroup& param)
|
|
{
|
|
return StringUtils::VectorValue::
|
|
get<std::string>(param.get<std::string>("quant"));
|
|
}
|
|
|
|
int numDigits(const std::vector<int>& steps)
|
|
{
|
|
if (steps.empty()) {
|
|
return 1;
|
|
}
|
|
|
|
const auto m =
|
|
*std::max_element(std::begin(steps), std::end(steps));
|
|
|
|
if (m == 0) {
|
|
return 1;
|
|
}
|
|
|
|
assert (m > 0);
|
|
|
|
return std::floor(std::log10(static_cast<double>(m))) + 1;
|
|
}
|
|
|
|
ReferenceSolution
|
|
loadReference(const ::Opm::ParameterGroup& param,
|
|
const std::string& quant,
|
|
const int step,
|
|
const int nDigits)
|
|
{
|
|
namespace fs = boost::filesystem;
|
|
|
|
using VRef = std::reference_wrapper<std::vector<double>>;
|
|
|
|
auto x = ReferenceSolution{};
|
|
auto ref = std::array<VRef,2>{{ std::ref(x.raw) ,
|
|
std::ref(x.SI ) }};
|
|
|
|
auto i = 0;
|
|
|
|
for (const auto* q : { "raw", "SI" }) {
|
|
auto fname = fs::path(param.get<std::string>("ref-dir"))
|
|
/ boost::algorithm::to_lower_copy(quant);
|
|
{
|
|
std::ostringstream os;
|
|
|
|
os << q << '-'
|
|
<< std::setw(nDigits) << std::setfill('0')
|
|
<< step << ".txt";
|
|
|
|
fname /= os.str();
|
|
}
|
|
|
|
fs::ifstream input(fname);
|
|
|
|
if (input) {
|
|
ref[i].get().assign(std::istream_iterator<double>(input),
|
|
std::istream_iterator<double>());
|
|
}
|
|
|
|
i += 1;
|
|
}
|
|
|
|
if (x.raw.size() != x.SI.size()) {
|
|
std::ostringstream os;
|
|
|
|
os << "Unable to Read Consistent Reference Data From '"
|
|
<< param.get<std::string>("ref-dir") << "' In Step "
|
|
<< step;
|
|
|
|
throw std::out_of_range(os.str());
|
|
}
|
|
|
|
return x;
|
|
}
|
|
|
|
void computeErrors(const Reference& ref,
|
|
const Calculated& calc,
|
|
AggregateErrors& E)
|
|
{
|
|
const auto diff =
|
|
VectorDifference(calc.data, ref.data); // calc - ref
|
|
|
|
using Vector1 = std::decay<decltype(diff)>::type;
|
|
using Vector2 = std::decay<decltype(ref.data)>::type;
|
|
using Ratio = VectorRatio<Vector1, Vector2>;
|
|
|
|
const auto rat = Ratio(diff, ref.data); // (tof - ref) / ref
|
|
|
|
auto abs = ErrorMeasurement{};
|
|
{
|
|
abs.volume = volumeMetric(diff);
|
|
abs.inf = pointMetric (diff);
|
|
}
|
|
|
|
auto rel = ErrorMeasurement{};
|
|
{
|
|
rel.volume = volumeMetric(rat);
|
|
rel.inf = pointMetric (rat);
|
|
}
|
|
|
|
E.absolute.push_back(std::move(abs));
|
|
E.relative.push_back(std::move(rel));
|
|
}
|
|
|
|
std::unique_ptr<Opm::ECLRestartData>
|
|
openRestartSet(const Opm::ECLCaseUtilities::ResultSet& rset,
|
|
const int step)
|
|
{
|
|
return std::unique_ptr<Opm::ECLRestartData> {
|
|
new Opm::ECLRestartData(rset.restartFile(step))
|
|
};
|
|
}
|
|
|
|
std::array<AggregateErrors, 2>
|
|
sampleDifferences(const ::Opm::ECLGraph& graph,
|
|
const ::Opm::ECLCaseUtilities::ResultSet& rset,
|
|
const ::Opm::ParameterGroup& param,
|
|
const std::string& quant,
|
|
const std::vector<int>& steps)
|
|
{
|
|
const auto ECLquant = boost::algorithm::to_upper_copy(quant);
|
|
|
|
auto unit = VectorUnits()
|
|
.getUnit(boost::algorithm::to_lower_copy(quant));
|
|
|
|
const auto start =
|
|
std::vector<Opm::FlowDiagnostics::CellSet>{};
|
|
|
|
const auto nDigits = numDigits(steps);
|
|
|
|
auto rstrt = std::unique_ptr<Opm::ECLRestartData>{};
|
|
|
|
auto E = std::array<AggregateErrors, 2>{};
|
|
|
|
for (const auto& step : steps) {
|
|
if (! (rset.isUnifiedRestart() && bool(rstrt))) {
|
|
// Separate (not unified) restart file or this is the first
|
|
// time we're selecting a report step.
|
|
rstrt = openRestartSet(rset, step);
|
|
}
|
|
|
|
if (! rstrt->selectReportStep(step)) {
|
|
continue;
|
|
}
|
|
|
|
const auto ref = loadReference(param, quant, step, nDigits);
|
|
|
|
{
|
|
const auto raw = Calculated {
|
|
graph.rawLinearisedCellData<double>(*rstrt, ECLquant)
|
|
};
|
|
|
|
computeErrors(Reference{ ref.raw }, raw, E[0]);
|
|
}
|
|
|
|
{
|
|
const auto SI = Calculated {
|
|
graph.linearisedCellData(*rstrt, ECLquant, unit)
|
|
};
|
|
|
|
computeErrors(Reference{ ref.SI }, SI, E[1]);
|
|
}
|
|
}
|
|
|
|
return E;
|
|
}
|
|
|
|
bool errorAcceptable(const std::vector<ErrorMeasurement>& E,
|
|
const double tol)
|
|
{
|
|
return std::accumulate(std::begin(E), std::end(E), true,
|
|
[tol](const bool ok, const ErrorMeasurement& e)
|
|
{
|
|
// Fine if at least one of .volume or .inf <= tol.
|
|
return ok && ! ((e.volume > tol) && (e.inf > tol));
|
|
});
|
|
}
|
|
|
|
bool everythingFine(const AggregateErrors& E,
|
|
const ErrorTolerance& tol)
|
|
{
|
|
return errorAcceptable(E.absolute, tol.absolute)
|
|
&& errorAcceptable(E.relative, tol.relative);
|
|
}
|
|
|
|
::Opm::ECLGraph
|
|
constructGraph(const Opm::ECLCaseUtilities::ResultSet& rset)
|
|
{
|
|
const auto I = ::Opm::ECLInitFileData(rset.initFile());
|
|
|
|
return ::Opm::ECLGraph::load(rset.gridFile(), I);
|
|
}
|
|
} // namespace Anonymous
|
|
|
|
int main(int argc, char* argv[])
|
|
try {
|
|
const auto prm = example::initParam(argc, argv);
|
|
const auto rset = example::identifyResultSet(prm);
|
|
const auto tol = testTolerances(prm);
|
|
|
|
const auto steps = rset.reportStepIDs();
|
|
const auto graph = constructGraph(rset);
|
|
|
|
auto all_ok = true;
|
|
for (const auto& quant : testQuantities(prm)) {
|
|
const auto E =
|
|
sampleDifferences(graph, rset, prm, quant, steps);
|
|
|
|
const auto ok =
|
|
everythingFine(E[0], tol) && everythingFine(E[1], tol);
|
|
|
|
std::cout << quant << ": " << (ok ? "OK" : "FAIL") << '\n';
|
|
|
|
all_ok = all_ok && ok;
|
|
}
|
|
|
|
if (! all_ok) {
|
|
return EXIT_FAILURE;
|
|
}
|
|
}
|
|
catch (const std::exception& e) {
|
|
std::cerr << "Caught Exception: " << e.what() << '\n';
|
|
|
|
return EXIT_FAILURE;
|
|
}
|