Merge pull request #638 from joakim-hove/udq-functions

Udq functions
This commit is contained in:
Joakim Hove 2019-02-25 09:35:00 +01:00 committed by GitHub
commit 40aad6d64e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 1143 additions and 12 deletions

View File

@ -131,6 +131,8 @@ if(ENABLE_ECL_INPUT)
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQEnums.cpp
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQInput.cpp
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQContext.cpp
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunction.cpp
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunctionTable.cpp
src/opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQExpression.cpp
src/opm/parser/eclipse/EclipseState/Schedule/VFPInjTable.cpp
src/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.cpp
@ -522,6 +524,8 @@ if(ENABLE_ECL_INPUT)
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQInput.hpp
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQWellSet.hpp
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQSet.hpp
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunction.hpp
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunctionTable.hpp
opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQExpression.hpp
opm/parser/eclipse/Deck/DeckItem.hpp
opm/parser/eclipse/Deck/Deck.hpp

View File

@ -0,0 +1,114 @@
/*
Copyright 2019 Equinor 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/>.
*/
#ifndef UDQFUNCTION_HPP
#define UDQFUNCTION_HPP
#include <stdexcept>
#include <vector>
#include <functional>
#include <random>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQSet.hpp>
namespace Opm {
class UDQFunction {
public:
explicit UDQFunction(const std::string& name);
virtual ~UDQFunction() = default;
const std::string& name() const;
private:
std::string m_name;
};
class UDQScalarFunction : public UDQFunction {
public:
UDQScalarFunction(const std::string&name, std::function<UDQScalar(const UDQSet& arg)> f);
UDQScalar eval(const UDQSet& arg) const;
static UDQScalar SUM(const UDQSet& arg);
static UDQScalar AVEA(const UDQSet& arg);
static UDQScalar AVEG(const UDQSet& arg);
static UDQScalar AVEH(const UDQSet& arg);
static UDQScalar MIN(const UDQSet& arg);
static UDQScalar MAX(const UDQSet& arg);
static UDQScalar NORM1(const UDQSet& arg);
static UDQScalar NORM2(const UDQSet& arg);
static UDQScalar NORMI(const UDQSet& arg);
static UDQScalar PROD(const UDQSet& arg);
private:
std::function<UDQScalar(const UDQSet& arg)> func;
};
class UDQUnaryElementalFunction : public UDQFunction {
public:
UDQUnaryElementalFunction(const std::string&name, std::function<UDQSet(const UDQSet& arg)> f);
UDQSet eval(const UDQSet& arg) const;
static UDQSet ABS(const UDQSet& arg);
static UDQSet DEF(const UDQSet& arg);
static UDQSet EXP(const UDQSet& arg);
static UDQSet IDV(const UDQSet& arg);
static UDQSet LN(const UDQSet& arg);
static UDQSet LOG(const UDQSet& arg);
static UDQSet NINT(const UDQSet& arg);
static UDQSet SORTA(const UDQSet& arg);
static UDQSet SORTD(const UDQSet& arg);
static UDQSet UNDEF(const UDQSet& arg);
static UDQSet RANDN(std::mt19937& rng, const UDQSet& arg);
static UDQSet RANDU(std::mt19937& rng, const UDQSet& arg);
private:
std::function<UDQSet(const UDQSet& arg)> func;
};
class UDQBinaryFunction : public UDQFunction {
public:
UDQBinaryFunction(const std::string& name, std::function<UDQSet(const UDQSet& lhs, const UDQSet& rhs)> f);
UDQSet eval(const UDQSet&, const UDQSet& arg) const;
static UDQSet EQ(double eps, const UDQSet& lhs, const UDQSet& rhs);
static UDQSet NE(double eps, const UDQSet& lhs, const UDQSet& rhs);
static UDQSet LE(double eps, const UDQSet& lhs, const UDQSet& rhs);
static UDQSet GE(double eps, const UDQSet& lhs, const UDQSet& rhs);
static UDQSet POW(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet LT(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet GT(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet ADD(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet MUL(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet SUB(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet DIV(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet UADD(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet UMUL(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet UMAX(const UDQSet& lhs, const UDQSet& rhs);
static UDQSet UMIN(const UDQSet& lhs, const UDQSet& rhs);
private:
std::function<UDQSet(const UDQSet& lhs, const UDQSet& rhs)> func;
};
}
#endif

View File

@ -0,0 +1,42 @@
/*
Copyright 2019 Equinor 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/>.
*/
#ifndef UDQFUNCTIONTABLE_HPP
#define UDQFUNCTIONTABLE_HPP
#include <unordered_map>
#include <memory>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunction.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQParams.hpp>
namespace Opm {
class UDQFunctionTable {
public:
UDQFunctionTable(UDQParams& params);
bool has_function(const std::string& name) const;
const UDQFunction& get(const std::string& name) const;
private:
void insert_function(std::unique_ptr<const UDQFunction> func);
UDQParams params;
std::unordered_map<std::string, std::unique_ptr<const UDQFunction>> function_table;
};
}
#endif

View File

@ -20,6 +20,7 @@
#ifndef OPM_UDQ_PARAMS_HPP
#define OPM_UDQ_PARAMS_HPP
#include <random>
namespace Opm {
@ -31,18 +32,22 @@ namespace Opm {
explicit UDQParams(const Deck& deck);
UDQParams();
bool reseedRNG() const noexcept;
int randomSeed() const noexcept;
void reseedRNG(int seed);
double range() const noexcept;
double undefinedValue() const noexcept;
double cmpEpsilon() const noexcept;
std::mt19937& sim_rng();
std::mt19937& true_rng();
private:
bool reseed_rng;
int random_seed;
double value_range;
double undefined_value;
double cmp_eps;
std::mt19937 m_sim_rng; // The sim_rng is seeded deterministiaclly at simulation start.
std::mt19937 m_true_rng; // The true_rng is seeded with a "true" random seed; this rng can be reset with reseedRNG()
};
}

View File

@ -0,0 +1,499 @@
/*
Copyright 2019 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 <unordered_set>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <random>
#include <numeric>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunction.hpp>
namespace Opm {
UDQFunction::UDQFunction(const std::string& name) :
m_name(name)
{
}
const std::string& UDQFunction::name() const {
return this->m_name;
}
UDQScalarFunction::UDQScalarFunction(const std::string&name, std::function<UDQScalar(const UDQSet& arg)> f) :
UDQFunction(name),
func(std::move(f))
{
}
UDQScalar UDQScalarFunction::eval(const UDQSet& arg) const {
return this->func(arg);
}
UDQSet UDQUnaryElementalFunction::eval(const UDQSet& arg) const {
return this->func(arg);
}
/*****************************************************************/
UDQScalar UDQScalarFunction::MIN(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
if (defined_values.empty())
return {};
return *std::min_element(defined_values.begin(), defined_values.end());
}
UDQScalar UDQScalarFunction::MAX(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
if (defined_values.empty())
return {};
return *std::max_element(defined_values.begin(), defined_values.end());
}
UDQScalar UDQScalarFunction::SUM(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
return std::accumulate(defined_values.begin(), defined_values.end(), 0.0);
}
UDQScalar UDQScalarFunction::PROD(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
return std::accumulate(defined_values.begin(), defined_values.end(), 1.0, std::multiplies<double>{});
}
UDQScalar UDQScalarFunction::AVEA(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
if (defined_values.empty())
return {};
return std::accumulate( defined_values.begin(), defined_values.end(), 0.0) / defined_values.size();
}
UDQScalar UDQScalarFunction::AVEG(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
if (defined_values.empty())
return {};
if (std::find_if(defined_values.begin(), defined_values.end(), [](double x) { return x <= 0; }) != defined_values.end())
throw std::invalid_argument("Function AVEG must have only positive arguments");
double log_mean = std::accumulate( defined_values.begin(), defined_values.end(), 0.0, [](double x, double y) { return x + std::log(std::fabs(y)); }) / defined_values.size();
return std::exp( log_mean );
}
UDQScalar UDQScalarFunction::AVEH(const UDQSet& arg) {
std::vector<double> defined_values = arg.defined_values();
if (defined_values.empty())
return {};
return defined_values.size() / std::accumulate(defined_values.begin(), defined_values.end(), 0.0, [](double x, double y) { return x + 1.0/y; });
}
UDQScalar UDQScalarFunction::NORM1(const UDQSet& arg) {
auto defined_values = arg.defined_values();
return std::accumulate( defined_values.begin(), defined_values.end(), 0.0, [](double x, double y) { return x + std::fabs(y);});
}
UDQScalar UDQScalarFunction::NORM2(const UDQSet& arg) {
auto defined_values = arg.defined_values();
return std::sqrt(std::inner_product(defined_values.begin(), defined_values.end(), defined_values.begin(), 0.0));
}
UDQScalar UDQScalarFunction::NORMI(const UDQSet& arg) {
auto defined_values = arg.defined_values();
return std::accumulate( defined_values.begin(), defined_values.end(), 0.0, [](const double x, const double y) { return std::max(x, std::fabs(y));});
}
UDQUnaryElementalFunction::UDQUnaryElementalFunction(const std::string&name, std::function<UDQSet(const UDQSet& arg)> f) :
UDQFunction(name),
func(f)
{}
UDQSet UDQUnaryElementalFunction::ABS(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, std::fabs(udq_value.value()));
}
return result;
}
UDQSet UDQUnaryElementalFunction::DEF(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, 1 );
}
return result;
}
UDQSet UDQUnaryElementalFunction::UNDEF(const UDQSet& arg) {
UDQSet result(arg.size());
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = arg[index];
if (!udq_value)
result.assign( index, 1 );
}
return result;
}
UDQSet UDQUnaryElementalFunction::IDV(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, 1 );
else
result.assign(index , 0);
}
return result;
}
UDQSet UDQUnaryElementalFunction::EXP(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, std::exp(udq_value.value()) );
}
return result;
}
UDQSet UDQUnaryElementalFunction::NINT(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, std::nearbyint(udq_value.value()) );
}
return result;
}
UDQSet UDQUnaryElementalFunction::RANDN(std::mt19937& rng, const UDQSet& arg) {
auto result = arg;
std::normal_distribution<double> dist(0,1);
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, dist(rng) );
}
return result;
}
UDQSet UDQUnaryElementalFunction::RANDU(std::mt19937& rng, const UDQSet& arg) {
auto result = arg;
std::uniform_real_distribution<double> dist(-1,1);
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value)
result.assign( index, dist(rng) );
}
return result;
}
UDQSet UDQUnaryElementalFunction::LN(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value) {
double elm = udq_value.value();
if (elm > 0)
result.assign(index, std::log(elm));
else
throw std::invalid_argument("Argument: " + std::to_string(elm) + " invalid for function LN");
}
}
return result;
}
UDQSet UDQUnaryElementalFunction::LOG(const UDQSet& arg) {
auto result = arg;
for (std::size_t index=0; index < result.size(); index++) {
auto& udq_value = result[index];
if (udq_value) {
double elm = udq_value.value();
if (elm > 0)
result.assign(index, std::log10(elm));
else
throw std::invalid_argument("Argument: " + std::to_string(elm) + " invalid for function LOG");
}
}
return result;
}
namespace {
UDQSet udq_union(const UDQSet& arg1, const UDQSet& arg2) {
if (arg1.size() != arg2.size())
throw std::invalid_argument("UDQ sets have incoimpatible size");
UDQSet result(arg1.size());
for (std::size_t index = 0; index < result.size(); index++) {
const auto& elm1 = arg1[index];
const auto& elm2 = arg2[index];
if (elm1.defined() != elm2.defined()) {
if (elm1)
result.assign(index, elm1.value());
if (elm2)
result.assign(index, elm2.value());
}
}
return result;
}
UDQSet udq_sort(const UDQSet& arg, std::size_t defined_size, const std::function<bool(int,int)>& cmp) {
auto result = arg;
std::vector<int> index(defined_size);
std::iota(index.begin(), index.end(), 0);
std::sort(index.begin(), index.end(), cmp);
std::size_t output_index = 0;
for (auto sort_index : index) {
while (!result[output_index])
output_index++;
result.assign(output_index, sort_index);
output_index++;
}
return result;
}
}
UDQSet UDQUnaryElementalFunction::SORTA(const UDQSet& arg) {
auto defined_values = arg.defined_values();
return udq_sort(arg, defined_values.size(), [&defined_values](int a, int b){ return defined_values[a] < defined_values[b]; });
}
UDQSet UDQUnaryElementalFunction::SORTD(const UDQSet& arg) {
auto defined_values = arg.defined_values();
return udq_sort(arg, defined_values.size(), [&defined_values](int a, int b){ return defined_values[a] > defined_values[b]; });
}
UDQBinaryFunction::UDQBinaryFunction(const std::string& name, std::function<UDQSet(const UDQSet& lhs, const UDQSet& rhs)> f) :
UDQFunction(name),
func(std::move(f))
{
}
UDQSet UDQBinaryFunction::eval(const UDQSet& lhs, const UDQSet& rhs) const {
return this->func(lhs, rhs);
}
UDQSet UDQBinaryFunction::LE(double eps, const UDQSet& lhs, const UDQSet& rhs) {
auto result = lhs - rhs;
auto rel_diff = result / lhs;
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm) {
double diff = rel_diff[index].value();
if (diff <= eps)
result.assign(index, 1);
else
result.assign(index, 0);
}
}
return result;
}
UDQSet UDQBinaryFunction::GE(double eps, const UDQSet& lhs, const UDQSet& rhs) {
auto result = lhs - rhs;
auto rel_diff = result / lhs;
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm) {
double diff = rel_diff[index].value();
if (diff >= -eps)
result.assign(index, 1);
else
result.assign(index, 0);
}
}
return result;
}
UDQSet UDQBinaryFunction::EQ(double eps, const UDQSet& lhs, const UDQSet& rhs) {
auto result = lhs - rhs;
auto rel_diff = result / lhs;
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm) {
double diff = std::fabs(rel_diff[index].value());
if (diff <= eps)
result.assign(index, 1);
else
result.assign(index, 0);
}
}
return result;
}
UDQSet UDQBinaryFunction::NE(double eps, const UDQSet& lhs, const UDQSet& rhs) {
auto result = UDQBinaryFunction::EQ(eps, lhs, rhs);
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm)
result.assign(index, 1 - elm.value());
}
return result;
}
UDQSet UDQBinaryFunction::GT(const UDQSet& lhs, const UDQSet& rhs) {
auto result = lhs - rhs;
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm) {
double diff = elm.value();
if (diff > 0)
result.assign(index, 1);
else
result.assign(index, 0);
}
}
return result;
}
UDQSet UDQBinaryFunction::LT(const UDQSet& lhs, const UDQSet& rhs) {
auto result = lhs - rhs;
for (std::size_t index=0; index < result.size(); index++) {
auto elm = result[index];
if (elm) {
double diff = elm.value();
if (diff < 0)
result.assign(index, 1);
else
result.assign(index, 0);
}
}
return result;
}
UDQSet UDQBinaryFunction::ADD(const UDQSet& lhs, const UDQSet& rhs) {
return lhs + rhs;
}
UDQSet UDQBinaryFunction::UADD(const UDQSet& lhs, const UDQSet& rhs) {
UDQSet result = udq_union(lhs,rhs);
for (std::size_t index=0; index < lhs.size(); index++) {
const auto& lhs_elm = lhs[index];
const auto& rhs_elm = rhs[index];
if (lhs_elm && rhs_elm)
result.assign(index, rhs_elm.value() + lhs_elm.value());
}
return result;
}
UDQSet UDQBinaryFunction::UMUL(const UDQSet& lhs, const UDQSet& rhs) {
UDQSet result = udq_union(lhs, rhs);
for (std::size_t index=0; index < lhs.size(); index++) {
const auto& lhs_elm = lhs[index];
const auto& rhs_elm = rhs[index];
if (lhs_elm && rhs_elm)
result.assign(index, rhs_elm.value() * lhs_elm.value());
}
return result;
}
UDQSet UDQBinaryFunction::UMIN(const UDQSet& lhs, const UDQSet& rhs) {
UDQSet result = udq_union(lhs, rhs);
for (std::size_t index=0; index < lhs.size(); index++) {
const auto& lhs_elm = lhs[index];
const auto& rhs_elm = rhs[index];
if (lhs_elm && rhs_elm)
result.assign(index, std::min(rhs_elm.value(), lhs_elm.value()));
}
return result;
}
UDQSet UDQBinaryFunction::UMAX(const UDQSet& lhs, const UDQSet& rhs) {
UDQSet result = udq_union(lhs, rhs);
for (std::size_t index=0; index < lhs.size(); index++) {
const auto& lhs_elm = lhs[index];
const auto& rhs_elm = rhs[index];
if (lhs_elm && rhs_elm)
result.assign(index, std::max(rhs_elm.value(), lhs_elm.value()));
}
return result;
}
UDQSet UDQBinaryFunction::MUL(const UDQSet& lhs, const UDQSet& rhs) {
return lhs * rhs;
}
UDQSet UDQBinaryFunction::SUB(const UDQSet& lhs, const UDQSet& rhs) {
return lhs - rhs;
}
UDQSet UDQBinaryFunction::DIV(const UDQSet& lhs, const UDQSet& rhs) {
return lhs / rhs;
}
UDQSet UDQBinaryFunction::POW(const UDQSet& lhs, const UDQSet& rhs) {
UDQSet result(lhs.size());
for (std::size_t index = 0; index < result.size(); index++) {
auto& lhs_elm = lhs[index];
auto& rhs_elm = rhs[index];
if (lhs_elm && rhs_elm)
result.assign(index, std::pow(lhs_elm.value(), rhs_elm.value()));
}
return result;
}
}

View File

@ -0,0 +1,110 @@
/*
Copyright 2019 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 <unordered_set>
#include <cmath>
#include <algorithm>
#include <functional>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunctionTable.hpp>
namespace Opm {
UDQFunctionTable::UDQFunctionTable(UDQParams& params) :
params(params)
{
// SCalar functions
this->insert_function( std::make_unique<UDQScalarFunction>("SUM", UDQScalarFunction::SUM) );
this->insert_function( std::make_unique<UDQScalarFunction>("AVEA", UDQScalarFunction::AVEA) );
this->insert_function( std::make_unique<UDQScalarFunction>("AVEG", UDQScalarFunction::AVEG) );
this->insert_function( std::make_unique<UDQScalarFunction>("AVEH", UDQScalarFunction::AVEH) );
this->insert_function( std::make_unique<UDQScalarFunction>("MAX", UDQScalarFunction::MAX) );
this->insert_function( std::make_unique<UDQScalarFunction>("MIN", UDQScalarFunction::MIN) );
this->insert_function( std::make_unique<UDQScalarFunction>("NORM1", UDQScalarFunction::NORM1) );
this->insert_function( std::make_unique<UDQScalarFunction>("NORM2", UDQScalarFunction::NORM2) );
this->insert_function( std::make_unique<UDQScalarFunction>("NORMI", UDQScalarFunction::NORMI) );
this->insert_function( std::make_unique<UDQScalarFunction>("PROD", UDQScalarFunction::PROD) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("ABS", UDQUnaryElementalFunction::ABS) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("DEF", UDQUnaryElementalFunction::DEF) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("EXP", UDQUnaryElementalFunction::EXP) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("IDV", UDQUnaryElementalFunction::IDV) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("LN", UDQUnaryElementalFunction::LN) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("LOG", UDQUnaryElementalFunction::LOG) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("NINT", UDQUnaryElementalFunction::NINT) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("SORTA", UDQUnaryElementalFunction::SORTA) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("SORTD", UDQUnaryElementalFunction::SORTD) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("UNDEF", UDQUnaryElementalFunction::UNDEF) );
const auto& randn = [ &sim_rng = this->params.sim_rng() ](const UDQSet& arg) { return UDQUnaryElementalFunction::RANDN(sim_rng, arg); };
const auto& randu = [ &sim_rng = this->params.sim_rng()] (const UDQSet& arg) { return UDQUnaryElementalFunction::RANDU(sim_rng, arg); };
const auto& true_randn = [ &true_rng = this->params.true_rng() ](const UDQSet& arg) { return UDQUnaryElementalFunction::RANDN(true_rng, arg); };
const auto& true_randu = [ &true_rng = this->params.true_rng() ](const UDQSet& arg) { return UDQUnaryElementalFunction::RANDU(true_rng, arg); };
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("RANDN", randn) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("RANDU", randu) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("RRNDN", true_randn) );
this->insert_function( std::make_unique<UDQUnaryElementalFunction>("RRNDU", true_randu) );
const auto& eq = [ eps = this->params.cmpEpsilon()](const UDQSet&lhs, const UDQSet&rhs) { return UDQBinaryFunction::EQ(eps, lhs, rhs); };
const auto& ne = [ eps = this->params.cmpEpsilon()](const UDQSet&lhs, const UDQSet&rhs) { return UDQBinaryFunction::NE(eps, lhs, rhs); };
const auto& ge = [ eps = this->params.cmpEpsilon()](const UDQSet&lhs, const UDQSet&rhs) { return UDQBinaryFunction::GE(eps, lhs, rhs); };
const auto& le = [ eps = this->params.cmpEpsilon()](const UDQSet&lhs, const UDQSet&rhs) { return UDQBinaryFunction::LE(eps, lhs, rhs); };
this->insert_function( std::make_unique<UDQBinaryFunction>("==", eq) );
this->insert_function( std::make_unique<UDQBinaryFunction>("!=", ne) );
this->insert_function( std::make_unique<UDQBinaryFunction>(">=", ge) );
this->insert_function( std::make_unique<UDQBinaryFunction>("<=", le) );
this->insert_function( std::make_unique<UDQBinaryFunction>("^", UDQBinaryFunction::POW ));
this->insert_function( std::make_unique<UDQBinaryFunction>("^", UDQBinaryFunction::POW ));
this->insert_function( std::make_unique<UDQBinaryFunction>("<", UDQBinaryFunction::LT ));
this->insert_function( std::make_unique<UDQBinaryFunction>(">", UDQBinaryFunction::GT ));
this->insert_function( std::make_unique<UDQBinaryFunction>("+", UDQBinaryFunction::ADD ));
this->insert_function( std::make_unique<UDQBinaryFunction>("*", UDQBinaryFunction::MUL ));
this->insert_function( std::make_unique<UDQBinaryFunction>("/", UDQBinaryFunction::DIV ));
this->insert_function( std::make_unique<UDQBinaryFunction>("-", UDQBinaryFunction::SUB ));
this->insert_function( std::make_unique<UDQBinaryFunction>("UADD", UDQBinaryFunction::UADD ));
this->insert_function( std::make_unique<UDQBinaryFunction>("UMUL", UDQBinaryFunction::UMUL ));
this->insert_function( std::make_unique<UDQBinaryFunction>("UMIN", UDQBinaryFunction::UMIN ));
this->insert_function( std::make_unique<UDQBinaryFunction>("UMAX", UDQBinaryFunction::UMAX ));
}
void UDQFunctionTable::insert_function(std::unique_ptr<const UDQFunction> func) {
auto name = func->name();
this->function_table.emplace( std::move(name), std::move(func) );
}
bool UDQFunctionTable::has_function(const std::string& name) const {
return this->function_table.count(name) > 0;
}
const UDQFunction& UDQFunctionTable::get(const std::string& name) const {
if (!this->has_function(name))
throw std::invalid_argument("No such function registered: " + name);
const auto& pair_ptr = this->function_table.find(name);
return *pair_ptr->second;
}
}

View File

@ -15,6 +15,7 @@
You should have received a copy of the GNU General Public License along with
OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <chrono>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/U.hpp>
@ -33,7 +34,25 @@ namespace Opm {
value_range(ParserKeywords::UDQPARAM::RANGE::defaultValue),
undefined_value(ParserKeywords::UDQPARAM::UNDEFINED_VALUE::defaultValue),
cmp_eps(ParserKeywords::UDQPARAM::CMP_EPSILON::defaultValue)
{}
{
/*
The std::random_device invokes instructions which are not properly
implemented on older valgrind versions, and this code will raise
SIGILL when running under valgrind:
https://stackoverflow.com/questions/37032339/valgrind-unrecognised-instruction
std::random_device r;
std::seed_seq seq{ r(), r(), r(), r() };
this->m_true_rng.seed( seq );
The random number stream is therefor seeded using time.
*/
auto now = std::chrono::high_resolution_clock::now();
auto ns = std::chrono::duration_cast<std::chrono::nanoseconds>(now.time_since_epoch());
this->m_true_rng.seed( ns.count() );
this->m_sim_rng.seed( this->random_seed );
}
UDQParams::UDQParams(const Deck& deck) :
@ -54,16 +73,21 @@ namespace Opm {
undefined_value = record.getItem("UNDEFINED_VALUE").get<double>(0);
cmp_eps = record.getItem("CMP_EPSILON").get<double>(0);
}
this->m_sim_rng.seed( this->random_seed );
}
bool UDQParams::reseedRNG() const noexcept {
return this->reseed_rng;
/*
If the internal flag reseed_rng is set to true, this method will reset the
rng true_rng with the seed applied; the purpose of this function is to be
able to ensure(?) that a restarted run gets the same sequence of random
numbers as the original run.
*/
void UDQParams::reseedRNG(int seed) {
if (this->reseed_rng)
this->m_true_rng.seed( seed );
}
int UDQParams::randomSeed() const noexcept {
return this->random_seed;
}
double UDQParams::range() const noexcept {
return this->value_range;
@ -76,4 +100,13 @@ namespace Opm {
double UDQParams::cmpEpsilon() const noexcept {
return this->cmp_eps;
}
std::mt19937& UDQParams::sim_rng() {
return this->m_sim_rng;
}
std::mt19937& UDQParams::true_rng() {
return this->m_true_rng;
}
}

View File

@ -31,6 +31,8 @@
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQExpression.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQContext.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQAssign.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunction.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQFunctionTable.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/SummaryState.hpp>
using namespace Opm;
@ -53,7 +55,7 @@ BOOST_AUTO_TEST_CASE(KEYWORDS) {
RUNSPEC
UDQDIMS
10* 'Y'/
10* 'N'/
UDQPARAM
3* 0.25 /
@ -65,8 +67,17 @@ UDQPARAM
auto runspec = Runspec(deck);
auto udq_params = runspec.udqParams();
BOOST_CHECK(udq_params.reseedRNG());
BOOST_CHECK_EQUAL(0.25, udq_params.cmpEpsilon());
// The reseed parameter is set to false, so the repeated calls to reseedRNG() should have
// no effect.
udq_params.reseedRNG(100);
auto r1 = udq_params.true_rng()();
udq_params.reseedRNG(100);
auto r2 = udq_params.true_rng()();
BOOST_CHECK( r1 != r2 );
}
@ -123,8 +134,21 @@ UDQ
Parser parser;
auto deck = parser.parseString(input);
auto udq_params = UDQParams(deck);
BOOST_CHECK_EQUAL(0.25, udq_params.cmpEpsilon());
auto udq_params1 = UDQParams(deck);
BOOST_CHECK_EQUAL(0.25, udq_params1.cmpEpsilon());
auto& sim_rng1 = udq_params1.sim_rng();
auto& true_rng1 = udq_params1.true_rng();
auto udq_params2 = UDQParams(deck);
auto& sim_rng2 = udq_params2.sim_rng();
auto& true_rng2 = udq_params2.true_rng();
BOOST_CHECK( sim_rng1() == sim_rng2() );
BOOST_CHECK( true_rng1() != true_rng2() );
udq_params1.reseedRNG(100);
udq_params2.reseedRNG(100);
BOOST_CHECK( true_rng1() == true_rng2() );
}
BOOST_AUTO_TEST_CASE(UDQ_CHANGE_UNITS_ILLEGAL) {
@ -315,6 +339,306 @@ BOOST_AUTO_TEST_CASE(UDQ_SET) {
}
BOOST_AUTO_TEST_CASE(UDQ_FUNCTION_TABLE) {
UDQParams params;
UDQFunctionTable udqft(params);
BOOST_CHECK(udqft.has_function("SUM"));
BOOST_CHECK(!udqft.has_function("NO_SUCH_FUNCTION"));
UDQSet arg(5);
arg.assign(0,1);
arg.assign(2,2);
arg.assign(4,4);
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("SUM"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 7);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("NORM1"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 7);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("NORM2"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), std::sqrt(1 + 4+ 16));
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("NORMI"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 4);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("MIN"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 1);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("MAX"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 4);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("AVEA"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 7.0/3);
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("AVEG"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), std::exp((std::log(1) + std::log(2.0) + std::log(4))/3));
}
{
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("PROD"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL(result.value(), 8.0);
}
{
UDQSet arg2(4);
arg2.assign(0,1);
arg2.assign(2,4);
arg2.assign(3,4);
const auto& func = dynamic_cast<const UDQScalarFunction&>(udqft.get("AVEH"));
auto result = func.eval(arg2);
BOOST_CHECK_EQUAL(result.value(), 2.0);
}
}
BOOST_AUTO_TEST_CASE(CMP_FUNCTIONS) {
UDQParams params;
UDQFunctionTable udqft(params);
UDQSet arg1(5);
UDQSet arg2(5);
UDQSet arg3(3);
arg1.assign(1,1);
arg1.assign(0,1);
arg1.assign(2,2);
arg1.assign(4,4);
arg2.assign(0, 0.9);
arg2.assign(2, 2.5);
arg2.assign(4, 4.0);
BOOST_CHECK_THROW(UDQBinaryFunction::EQ(0.25, arg1, arg3), std::invalid_argument);
{
auto result = UDQBinaryFunction::EQ(0, arg1, arg2);
BOOST_CHECK_EQUAL( result.defined_size(), 3 );
BOOST_CHECK_EQUAL( result[0].value(), 0);
BOOST_CHECK_EQUAL( result[2].value(), 0);
BOOST_CHECK_EQUAL( result[4].value(), 1);
result = UDQBinaryFunction::EQ(0.20, arg1, arg2);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK_EQUAL( result[2].value(), 0);
BOOST_CHECK_EQUAL( result[4].value(), 1);
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get("=="));
result = func.eval(arg1, arg2);
BOOST_CHECK_EQUAL( result[0].value(), 0);
BOOST_CHECK_EQUAL( result[2].value(), 0);
BOOST_CHECK_EQUAL( result[4].value(), 1);
}
{
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get("<"));
auto result = func.eval(arg1, arg2);
BOOST_CHECK_EQUAL( result.defined_size(), 3 );
BOOST_CHECK_EQUAL( result[0].value(), 0);
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK_EQUAL( result[4].value(), 0);
}
{
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get(">"));
auto result = func.eval(arg1, arg2);
BOOST_CHECK_EQUAL( result.defined_size(), 3 );
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK_EQUAL( result[2].value(), 0);
BOOST_CHECK_EQUAL( result[4].value(), 0);
}
{
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get("^"));
UDQSet arg1(4);
UDQSet arg2(4);
for (std::size_t i=0; i < arg1.size(); i++) {
arg1.assign(i, i + 1);
arg2.assign(i, 2);
}
auto result = func.eval(arg1, arg2);
for (std::size_t i=0; i < arg1.size(); i++)
BOOST_CHECK_EQUAL( result[i].value(), (i+1)*(i+1));
}
{
auto result = UDQBinaryFunction::GE(1.0, arg1, arg2);
BOOST_CHECK_EQUAL( result[0].value(), 1);
// This is bisarre - but due to the large epsilon 2 and 2.5 compare as
// equal; and then we evaluate 2 >= 2.5 as TRUE!
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK_EQUAL( result[4].value(), 1);
}
{
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get("<="));
auto result = func.eval(arg1, arg2);
BOOST_CHECK_EQUAL( result[0].value(), 0);
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK_EQUAL( result[4].value(), 1);
}
}
BOOST_AUTO_TEST_CASE(BAD_CAST) {
UDQParams params;
UDQFunctionTable udqft(params);
BOOST_CHECK_THROW( dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("==")), std::bad_cast);
}
BOOST_AUTO_TEST_CASE(ELEMENTAL_UNARY_FUNCTIONS) {
UDQParams params;
UDQFunctionTable udqft(params);
UDQSet arg(5);
arg.assign(0,1);
arg.assign(2,2);
arg.assign(4,4);
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("ABS"));
UDQSet arg2(5);
arg2.assign(0,1);
arg2.assign(2,-2);
arg2.assign(4,4);
auto result = func.eval(arg2);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK_EQUAL( result[2].value(), 2);
BOOST_CHECK_EQUAL( result[4].value(), 4);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("DEF"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK_EQUAL( result[4].value(), 1);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("UNDEF"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[1].value(), 1);
BOOST_CHECK_EQUAL( result[3].value(), 1);
BOOST_CHECK_EQUAL( result.defined_size(), 2);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("EXP"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), std::exp(1));
BOOST_CHECK_EQUAL( result[2].value(), std::exp(2));
BOOST_CHECK_EQUAL( result[4].value(), std::exp(4));
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("IDV"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK_EQUAL( result[1].value(), 0);
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK_EQUAL( result[3].value(), 0);
BOOST_CHECK_EQUAL( result[4].value(), 1);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("LOG"));
UDQSet arg(3);
arg.assign(0, 10);
arg.assign(2,1000);
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK( !result[1] );
BOOST_CHECK_EQUAL( result[2].value(), 3);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("NINT"));
UDQSet arg(3);
arg.assign(0, 0.75);
arg.assign(2, 1.25);
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 1);
BOOST_CHECK( !result[1] );
BOOST_CHECK_EQUAL( result[2].value(), 1);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("RANDN"));
UDQSet arg(3);
arg.assign(0, -1.0);
arg.assign(2, -1.0);
auto result1 = func.eval(arg);
auto result2 = func.eval(arg);
BOOST_CHECK( result1[0].value() != -1.0);
BOOST_CHECK( !result1[1] );
BOOST_CHECK( result1[2].value() != -1.0);
BOOST_CHECK( result1[0].value() != result2[0].value());
BOOST_CHECK( result1[2].value() != result2[2].value());
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("SORTA"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 0);
BOOST_CHECK( !result[1] );
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK( !result[3] );
BOOST_CHECK_EQUAL( result[4].value(), 2);
}
{
const auto& func = dynamic_cast<const UDQUnaryElementalFunction&>(udqft.get("SORTD"));
auto result = func.eval(arg);
BOOST_CHECK_EQUAL( result[0].value(), 2);
BOOST_CHECK( !result[1] );
BOOST_CHECK_EQUAL( result[2].value(), 1);
BOOST_CHECK( !result[3] );
BOOST_CHECK_EQUAL( result[4].value(), 0);
}
}
BOOST_AUTO_TEST_CASE(UNION_FUNCTIONS) {
UDQParams params;
UDQFunctionTable udqft(params);
UDQSet arg1(5);
UDQSet arg2(5);
arg1.assign(0,1);
arg1.assign(2,2);
arg2.assign(0, 1.0);
arg2.assign(3, 3 );
const auto& func = dynamic_cast<const UDQBinaryFunction&>(udqft.get("UADD"));
auto result = func.eval(arg1, arg2);
BOOST_CHECK_EQUAL( 3, result.defined_size() );
BOOST_CHECK_EQUAL( 2, result[0].value() );
BOOST_CHECK_EQUAL( 2, result[2].value() );
BOOST_CHECK_EQUAL( 3, result[3].value() );
}
BOOST_AUTO_TEST_CASE(FUNCTIONS_INVALID_ARGUMENT) {
UDQSet arg(3);
arg.assign(0, -1);
BOOST_REQUIRE_THROW( UDQScalarFunction::AVEG(arg), std::invalid_argument);
BOOST_REQUIRE_THROW( UDQUnaryElementalFunction::LOG(arg), std::invalid_argument);
BOOST_REQUIRE_THROW( UDQUnaryElementalFunction::LN(arg), std::invalid_argument);
}
BOOST_AUTO_TEST_CASE(UDQ_SET_DIV) {
UDQSet s(5);
s.assign(0,1);