opm-simulators/tests/test_region_phase_pvaverage.cpp
Bård Skaflestad d1d22668ea Add Utility for Region-Based Phase-Filled Averages
The initial use case is calculating the phase-filled pore-volume
weighted average of the fluid mass densities per PVT region.  This
value goes into calculating depth-corrected per-cell phase pressure
values such as the BPPO and BPPG summary vectors.

This class manages a single linear array which separately tracks the
averages' numerators and denominators as running sums per region and
region set.  We pick this data structure to simplify the cross-rank
reduction needed in MPI parallel runs.  Client code is expected to
add individual per-cell and per-phase contributions using the
addCell() member function and then call the accumulateParallel()
member to affect the cross-rank reduction.  The averages will then
be available through the fieldValue() and value() member functions.

As a further view towards the initial use case, we track two
different types of average per phase--one for the phase-filled
volume and one for the pore-volume filled volume.  The latter is the
average we would get for the case of the phase saturation being one
throughout the region.  This alternative value is the fallback
option for the case of the phase saturation being identically zero
throughout the region.
2024-02-28 12:08:24 +01:00

1069 lines
30 KiB
C++

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright 2024 Equinor.
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 <config.h>
#define BOOST_TEST_MODULE Region_Phase_PVAverage
#define BOOST_TEST_NO_MAIN
#include <boost/test/unit_test.hpp>
#include <opm/simulators/flow/RegionPhasePVAverage.hpp>
#include <opm/simulators/flow/partitionCells.hpp>
#include <opm/simulators/utils/ParallelCommunication.hpp>
#include <dune/common/parallel/mpihelper.hh>
#include <array>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iostream>
#include <iterator>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <unordered_map>
#include <vector>
#include <fmt/format.h>
namespace {
#if HAVE_MPI
struct MPIError
{
MPIError(std::string_view errstr, const int ec)
: errorstring { errstr }
, errorcode { ec }
{}
std::string errorstring;
int errorcode;
};
void MPI_err_handler(MPI_Comm*, int* err_code, ...)
{
std::array<char, MPI_MAX_ERROR_STRING> err_string_vec{'\0'};
auto err_length = 0;
MPI_Error_string(*err_code, err_string_vec.data(), &err_length);
auto err_string = std::string_view {
err_string_vec.data(),
static_cast<std::string_view::size_type>(err_length)
};
std::cerr << "An MPI Error ocurred:\n -> " << err_string << '\n';
throw MPIError { err_string, *err_code };
}
// Register a throwing error handler to allow for debugging with
//
// catch throw
//
// in GDB.
void register_error_handler()
{
MPI_Errhandler handler{};
MPI_Comm_create_errhandler(MPI_err_handler, &handler);
MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler);
}
#else // !HAVE_MPI
void register_error_handler()
{}
#endif // HAVE_MPI
template <typename GlobalVector>
auto selectSubset(const std::vector<bool>& keep,
const GlobalVector& global)
{
using VT = std::remove_cv_t<std::remove_reference_t<
typename std::iterator_traits<decltype(global.begin())>::value_type>>;
auto local = std::vector<VT>{};
const auto n = keep.size();
auto elm = global.begin();
for (auto i = 0*n; i < n; ++i, ++elm) {
if (keep[i]) {
local.push_back(*elm);
}
}
return local;
}
std::vector<bool> cellSubset(const int rank, const std::vector<int>& p)
{
auto keep = std::vector<bool>(p.size(), false);
const auto n = p.size();
for (auto i = 0*n; i < n; ++i) {
keep[i] = p[i] == rank;
}
return keep;
}
std::vector<int> fipnum()
{
return {
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
};
}
std::vector<int> fiplayer()
{
return {
0, 0, 0,
0, 0, 0,
0, 0, 0,
1, 1, 1,
1, 1, 1,
1, 1, 1,
2, 2, 2,
2, 2, 2,
2, 2, 2,
};
}
class RegionSets
{
public:
RegionSets();
RegionSets& selectSubset(const std::vector<bool>& keep);
std::vector<std::string> names() const;
const std::vector<int>& operator()(const std::string& rname) const;
private:
std::unordered_map<std::string, std::vector<int>> regions_;
};
RegionSets::RegionSets()
: regions_ {
{ "FIPNUM", fipnum() },
{ "FIPLRS", fiplayer() },
}
{}
RegionSets& RegionSets::selectSubset(const std::vector<bool>& keep)
{
for (auto& [name, region] : this->regions_) {
if (region.size() != keep.size()) {
continue;
}
auto local = ::selectSubset(keep, region);
region.swap(local);
}
return *this;
}
std::vector<std::string> RegionSets::names() const
{
auto rnames = std::vector<std::string> {};
rnames.reserve(this->regions_.size());
for (const auto& region : this->regions_) {
rnames.push_back(region.first);
}
return rnames;
}
const std::vector<int>& RegionSets::operator()(const std::string& rname) const
{
auto rpos = this->regions_.find(rname);
if (rpos == this->regions_.end()) {
throw std::invalid_argument {
fmt::format("Unknown region set '{}'", rname)
};
}
return rpos->second;
}
bool init_unit_test_func()
{
return true;
}
} // Anonymous namespace
// ===========================================================================
using PVAvg = Opm::RegionPhasePoreVolAverage;
BOOST_AUTO_TEST_SUITE(Sequential)
BOOST_AUTO_TEST_CASE(Single_Phase_Single_Val_Per_Layer)
{
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto rset = RegionSets{};
const auto numPhases = std::size_t{1};
const auto p = PVAvg::Phase {0};
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
BOOST_CHECK_MESSAGE(std::isnan(avgCalc.fieldValue(p)),
"Field value must be NaN prior to accumulation");
avgCalc.prepareAccumulation();
{
const auto activeCell = 0*9 + std::size_t{0};
const auto cv = PVAvg::CellValue {
1.0, 0.25, 200.0
};
avgCalc.addCell(activeCell, p, cv);
}
{
const auto activeCell = 1*9 + std::size_t{0};
const auto cv = PVAvg::CellValue {
2.0, 0.0, 200.0
};
avgCalc.addCell(activeCell, p, cv);
}
{
const auto activeCell = 2*9 + std::size_t{0};
const auto cv = PVAvg::CellValue {
3.0, 0.5, 200.0
};
avgCalc.addCell(activeCell, p, cv);
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(p), 7.0/3.0, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", p, PVAvg::Region{0}), 7.0/3.0, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{0}), 1.0, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{1}), 2.0, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{2}), 3.0, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Single_Phase_Full)
{
const auto x = std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
};
const auto s = std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
};
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto rset = RegionSets{};
const auto numPhases = std::size_t{1};
const auto p = PVAvg::Phase {0};
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], s[c], 200.0
};
avgCalc.addCell(c, p, cv);
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(p), 1.61, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", p, PVAvg::Region{0}), 1.61, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{0}), 1.11, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{1}), 2.11, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{2}), 3.11, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Single_Phase_Full_Varying_PV)
{
const auto x = std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
};
const auto s = std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
};
const auto pv = std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
};
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto rset = RegionSets{};
const auto numPhases = std::size_t{1};
const auto p = PVAvg::Phase {0};
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], s[c], pv[c]
};
avgCalc.addCell(c, p, cv);
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(p), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", p, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Three_Phase_All_Present_Full_Varying_PV)
{
const auto x = std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
};
const auto s = std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
};
const auto pv = std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
};
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto transform = std::array {
std::function { [](const double s0) { return s0 / 2; } },
std::function { [](const double s0) { return s0 / 4; } },
std::function { [](const double s0) { return 1 - (3*s0/4); } },
};
const auto rset = RegionSets{};
const auto numPhases = transform.size();
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
{
auto p = PVAvg::Phase {0};
for (const auto& t : transform) {
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], t(s[c]), pv[c]
};
avgCalc.addCell(c, p, cv);
}
++p.ix;
}
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{1}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{2}), 1.471079741628658, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{0}, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{1}, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{2}, PVAvg::Region{0}), 1.471079741628658, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{0}), 1.111326195193250, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{1}), 2.156805281711179, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Three_Phase_Two_Present_Full_Varying_PV)
{
const auto x = std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
};
const auto s = std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.25, 0.30, 0.25,
0.30, 0.35, 0.30,
0.25, 0.30, 0.25,
};
const auto pv = std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
};
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto transform = std::array {
std::function { [](const double s0) { return s0; } },
std::function { [](const double s0) { return 0 * s0; } },
std::function { [](const double s0) { return 1 - s0; } },
};
const auto rset = RegionSets{};
const auto numPhases = transform.size();
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
{
auto p = PVAvg::Phase {0};
for (const auto& t : transform) {
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], t(s[c]), pv[c]
};
avgCalc.addCell(c, p, cv);
}
++p.ix;
}
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{0}), 2.203870000146282, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{1}), 1.460875637211809, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{2}), 1.388846263879987, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{0}, PVAvg::Region{0}), 2.203870000146282, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{1}, PVAvg::Region{0}), 1.460875637211809, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{2}, PVAvg::Region{0}), 1.388846263879987, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{2}), 3.081133011812399, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{0}), 1.111895506705607, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{1}), 2.156118568232662, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{0}), 1.111126811726796, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{1}), 2.157055514795793, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{2}), 3.079851244928804, 1.0e-8);
}
BOOST_AUTO_TEST_SUITE_END() // Sequential
// ---------------------------------------------------------------------------
BOOST_AUTO_TEST_SUITE(Parallel)
BOOST_AUTO_TEST_CASE(Single_Phase_Full)
{
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto keep = cellSubset(comm.rank(), Opm::partitionCellsSimple(3 * 3 * 3, comm.size()).first);
const auto x = selectSubset(keep, std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
});
const auto s = selectSubset(keep, std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
});
const auto rset = RegionSets{}.selectSubset(keep);
const auto numPhases = std::size_t{1};
const auto p = PVAvg::Phase {0};
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], s[c], 200.0
};
avgCalc.addCell(c, p, cv);
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(p), 1.61, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", p, PVAvg::Region{0}), 1.61, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{0}), 1.11, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{1}), 2.11, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{2}), 3.11, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Single_Phase_Full_Varying_PV)
{
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto keep = cellSubset(comm.rank(), Opm::partitionCellsSimple(3 * 3 * 3, comm.size()).first);
const auto x = selectSubset(keep, std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
});
const auto s = selectSubset(keep, std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
});
const auto pv = selectSubset(keep, std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
});
const auto rset = RegionSets{}.selectSubset(keep);
const auto numPhases = std::size_t{1};
const auto p = PVAvg::Phase {0};
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], s[c], pv[c]
};
avgCalc.addCell(c, p, cv);
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(p), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", p, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", p, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Three_Phase_All_Present_Full_Varying_PV)
{
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto keep = cellSubset(comm.rank(), Opm::partitionCellsSimple(3 * 3 * 3, comm.size()).first);
const auto x = selectSubset(keep, std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
});
const auto s = selectSubset(keep, std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.0, 0.0, 0.0, // s=0 in layer 3 => use PV average instead
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
});
const auto pv = selectSubset(keep, std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
});
const auto transform = std::array {
std::function { [](const double s0) { return s0 / 2; } },
std::function { [](const double s0) { return s0 / 4; } },
std::function { [](const double s0) { return 1 - (3*s0/4); } },
};
const auto rset = RegionSets{}.selectSubset(keep);
const auto numPhases = transform.size();
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
{
auto p = PVAvg::Phase {0};
for (const auto& t : transform) {
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], t(s[c]), pv[c]
};
avgCalc.addCell(c, p, cv);
}
++p.ix;
}
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{1}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{2}), 1.471079741628658, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{0}, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{1}, PVAvg::Region{0}), 1.128401081038469, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{2}, PVAvg::Region{0}), 1.471079741628658, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{0}), 1.111326195193250, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{1}), 2.156805281711179, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
}
BOOST_AUTO_TEST_CASE(Three_Phase_Two_Present_Full_Varying_PV)
{
auto comm = Opm::Parallel::Communication {
Dune::MPIHelper::getCommunicator()
};
const auto keep = cellSubset(comm.rank(), Opm::partitionCellsSimple(3 * 3 * 3, comm.size()).first);
const auto x = selectSubset(keep, std::array {
// K=1
1.0, 1.01, 1.02,
1.1, 1.11, 1.12,
1.2, 1.21, 1.22,
// K=2
2.0, 2.01, 2.02,
2.1, 2.11, 2.12,
2.2, 2.21, 2.22,
// K=3
3.0, 3.01, 3.02,
3.1, 3.11, 3.12,
3.2, 3.21, 3.22,
});
const auto s = selectSubset(keep, std::array {
// K=1
0.0, 0.1, 0.0,
0.1, 0.4, 0.1,
0.0, 0.1, 0.0,
// K=2
0.1, 0.0, 0.1,
0.0, 0.4, 0.0,
0.1, 0.0, 0.1,
// K=3
0.25, 0.30, 0.25,
0.30, 0.35, 0.30,
0.25, 0.30, 0.25,
});
const auto pv = selectSubset(keep, std::array {
// K=1
1729.0, 271.82, 3141.5 ,
0.1, 1000.0, 321.09,
4321.0, 1234.5, 67.89,
// K=2
0.12, 0.34, 0.56,
0.78, 0.90, 1.23,
1.45, 1.67, 1.89,
// K=3
500.0, 450.0, 400.0,
350.0, 1.0, 300.0,
250.0, 200.0, 150.0,
});
const auto transform = std::array {
std::function { [](const double s0) { return s0; } },
std::function { [](const double s0) { return 0 * s0; } },
std::function { [](const double s0) { return 1 - s0; } },
};
const auto rset = RegionSets{}.selectSubset(keep);
const auto numPhases = transform.size();
auto avgCalc = PVAvg {
comm, numPhases, rset.names(), rset
};
avgCalc.prepareAccumulation();
{
auto p = PVAvg::Phase {0};
for (const auto& t : transform) {
for (auto nc = x.size(), c = 0*nc; c < nc; ++c) {
const auto cv = PVAvg::CellValue {
x[c], t(s[c]), pv[c]
};
avgCalc.addCell(c, p, cv);
}
++p.ix;
}
}
avgCalc.accumulateParallel();
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{0}), 2.203870000146282, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{1}), 1.460875637211809, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.fieldValue(PVAvg::Phase{2}), 1.388846263879987, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{0}, PVAvg::Region{0}), 2.203870000146282, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{1}, PVAvg::Region{0}), 1.460875637211809, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPNUM", PVAvg::Phase{2}, PVAvg::Region{0}), 1.388846263879987, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{0}), 1.127070395417597, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{1}), 2.146062992125984, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{0}, PVAvg::Region{2}), 3.081133011812399, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{0}), 1.111895506705607, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{1}), 2.156118568232662, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{1}, PVAvg::Region{2}), 3.080203767781622, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{0}), 1.111126811726796, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{1}), 2.157055514795793, 1.0e-8);
BOOST_CHECK_CLOSE(avgCalc.value("FIPLRS", PVAvg::Phase{2}, PVAvg::Region{2}), 3.079851244928804, 1.0e-8);
}
BOOST_AUTO_TEST_SUITE_END() // Parallel
// ===========================================================================
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
register_error_handler();
return boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv);
}