Merge pull request #3358 from bska/aqunum-initpress-array

Chase Type Specific Aquifer Data API Change
This commit is contained in:
Joakim Hove
2021-06-25 15:23:39 +02:00
committed by GitHub
8 changed files with 201 additions and 85 deletions

View File

@@ -37,6 +37,7 @@
#include <ebos/femcpgridcompat.hh>
#endif
#include <algorithm>
#include <cassert>
#include <stdexcept>
#include <string>
@@ -627,12 +628,20 @@ public:
this->unpackCommon(data, aq);
if (data.aquFet != nullptr) {
this->unpackAquFet(std::move(data.aquFet), aq);
if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
aquFet != nullptr)
{
this->unpackAquFet(*aquFet, aq);
}
if (data.aquCT != nullptr) {
this->unpackCarterTracy(std::move(data.aquCT), aq);
else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
aquCT != nullptr)
{
this->unpackCarterTracy(*aquCT, aq);
}
else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
aquNum != nullptr)
{
this->unpackNumericAquifer(*aquNum, aq);
}
}
}
@@ -648,14 +657,15 @@ private:
aq.volume += data.volume;
}
void unpackAquFet(std::shared_ptr<data::FetkovichData> aquFet, data::AquiferData& aq)
void unpackAquFet(const data::FetkovichData& aquFet, data::AquiferData& aq)
{
if (aq.aquFet == nullptr) {
aq.aquFet = std::move(aquFet);
if (! aq.typeData.is<data::AquiferType::Fetkovich>()) {
auto* tData = aq.typeData.create<data::AquiferType::Fetkovich>();
*tData = aquFet;
}
else {
const auto& src = *aquFet;
auto& dst = *aq.aquFet;
const auto& src = aquFet;
auto& dst = *aq.typeData.getMutable<data::AquiferType::Fetkovich>();
dst.initVolume = std::max(dst.initVolume , src.initVolume);
dst.prodIndex = std::max(dst.prodIndex , src.prodIndex);
@@ -663,19 +673,46 @@ private:
}
}
void unpackCarterTracy(std::shared_ptr<data::CarterTracyData> aquCT, data::AquiferData& aq)
void unpackCarterTracy(const data::CarterTracyData& aquCT, data::AquiferData& aq)
{
if (aq.aquCT == nullptr) {
aq.aquCT = std::move(aquCT);
if (! aq.typeData.is<data::AquiferType::CarterTracy>()) {
auto* tData = aq.typeData.create<data::AquiferType::CarterTracy>();
*tData = aquCT;
}
else {
const auto& src = *aquCT;
auto& dst = *aq.aquCT;
const auto& src = aquCT;
auto& dst = *aq.typeData.getMutable<data::AquiferType::CarterTracy>();
dst.timeConstant = std::max(dst.timeConstant , src.timeConstant);
dst.influxConstant = std::max(dst.influxConstant, src.influxConstant);
dst.waterDensity = std::max(dst.waterDensity , src.waterDensity);
dst.waterViscosity = std::max(dst.waterViscosity, src.waterViscosity);
dst.dimensionless_time = std::max(dst.dimensionless_time , src.dimensionless_time);
dst.dimensionless_pressure = std::max(dst.dimensionless_pressure, src.dimensionless_pressure);
}
}
void unpackNumericAquifer(const data::NumericAquiferData& aquNum, data::AquiferData& aq)
{
if (! aq.typeData.is<data::AquiferType::Numerical>()) {
auto* tData = aq.typeData.create<data::AquiferType::Numerical>();
*tData = aquNum;
}
else {
const auto& src = aquNum;
auto& dst = *aq.typeData.getMutable<data::AquiferType::Numerical>();
std::transform(src.initPressure.begin(),
src.initPressure.end(),
dst.initPressure.begin(),
dst.initPressure.begin(),
[](const double p0_1, const double p0_2)
{
return std::max(p0_1, p0_2);
});
}
}
};

View File

@@ -81,15 +81,14 @@ public:
}
data.volume = this->W_flux_.value();
data.initPressure = this->pa0_;
data.type = data::AquiferType::CarterTracy;
data.aquCT = std::make_shared<data::CarterTracyData>();
data.aquCT->timeConstant = this->aquct_data_.timeConstant();
data.aquCT->influxConstant = this->aquct_data_.influxConstant();
data.aquCT->waterDensity = this->aquct_data_.waterDensity();
data.aquCT->waterViscosity = this->aquct_data_.waterViscosity();
data.aquCT->dimensionless_time = this->dimensionless_time_;
data.aquCT->dimensionless_pressure = this->dimensionless_pressure_;
auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
aquCT->timeConstant = this->aquct_data_.timeConstant();
aquCT->influxConstant = this->aquct_data_.influxConstant();
aquCT->waterDensity = this->aquct_data_.waterDensity();
aquCT->waterViscosity = this->aquct_data_.waterViscosity();
aquCT->dimensionless_time = this->dimensionless_time_;
aquCT->dimensionless_pressure = this->dimensionless_pressure_;
return data;
}

View File

@@ -27,6 +27,7 @@ along with OPM. If not, see <http://www.gnu.org/licenses/>.
#include <exception>
#include <stdexcept>
#include <utility>
namespace Opm
{
@@ -82,12 +83,11 @@ public:
});
data.volume = this->W_flux_.value();
data.initPressure = this->pa0_;
data.type = data::AquiferType::Fetkovich;
data.aquFet = std::make_shared<data::FetkovichData>();
data.aquFet->initVolume = this->aqufetp_data_.initial_watvolume;
data.aquFet->prodIndex = this->aqufetp_data_.prod_index;
data.aquFet->timeConstant = this->aqufetp_data_.timeConstant();
auto* aquFet = data.typeData.template create<data::AquiferType::Fetkovich>();
aquFet->initVolume = this->aqufetp_data_.initial_watvolume;
aquFet->prodIndex = this->aqufetp_data_.prod_index;
aquFet->timeConstant = this->aqufetp_data_.timeConstant();
return data;
}
@@ -99,7 +99,7 @@ protected:
void assignRestartData(const data::AquiferData& xaq) override
{
if (xaq.type != data::AquiferType::Fetkovich) {
if (! xaq.typeData.is<data::AquiferType::Fetkovich>()) {
throw std::invalid_argument {
"Analytic aquifer data for unexpected aquifer "
"type passed to Fetkovich aquifer"

View File

@@ -24,6 +24,8 @@
#include <opm/output/data/Aquifer.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
#include <utility>
namespace Opm
{
template <typename TypeTag>
@@ -52,22 +54,22 @@ public:
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebos_simulator,
const int* global_cell)
: id_(aquifer.id())
, ebos_simulator_(ebos_simulator)
, flux_rate_(0.)
, cumulative_flux_(0.)
, global_cell_(global_cell)
: id_(aquifer.id())
, ebos_simulator_(ebos_simulator)
, flux_rate_(0.)
, cumulative_flux_(0.)
, global_cell_(global_cell)
, init_pressure_(aquifer.numCells(), 0.0)
{
this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < aquifer.numCells(); ++idx) {
const auto& cell = *(aquifer.getCellPrt(idx));
const int global_idx = cell.global_index;
const auto search = cartesian_to_compressed.find(global_idx);
const auto* cell = aquifer.getCellPrt(idx);
// Due to parallelisation, the cell might not exist in the current process
auto search = cartesian_to_compressed.find(cell->global_index);
if (search != cartesian_to_compressed.end()) {
const int cell_idx = cartesian_to_compressed.at(global_idx);
this->cell_to_aquifer_cell_idx_[cell_idx] = idx;
this->cell_to_aquifer_cell_idx_[search->second] = idx;
}
}
}
@@ -87,19 +89,20 @@ public:
data::AquiferData aquiferData() const
{
data::AquiferData data;
data.aquiferID = this->id_;
data.initPressure = this->init_pressure_;
data.aquiferID = this->aquiferID();
data.pressure = this->pressure_;
data.fluxRate = this->flux_rate_;
data.volume = this->cumulative_flux_;
data.type = data::AquiferType::Numerical;
auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
aquNum->initPressure = this->init_pressure_;
return data;
}
void initialSolutionApplied()
{
this->init_pressure_ = this->calculateAquiferPressure();
this->pressure_ = this->init_pressure_;
this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
this->flux_rate_ = 0.;
this->cumulative_flux_ = 0.;
}
@@ -115,13 +118,19 @@ private:
double flux_rate_; // aquifer influx rate
double cumulative_flux_; // cumulative aquifer influx
const int* global_cell_; // mapping to global index
double init_pressure_;
std::vector<double> init_pressure_{};
double pressure_; // aquifer pressure
// TODO: maybe unordered_map can also do the work to save memory?
std::vector<int> cell_to_aquifer_cell_idx_;
double calculateAquiferPressure() const
{
auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
return this->calculateAquiferPressure(capture);
}
double calculateAquiferPressure(std::vector<double>& cell_pressure) const
{
double sum_pressure_watervolume = 0.;
double sum_watervolume = 0.;
@@ -158,11 +167,17 @@ private:
const double water_volume = volume * porosity * water_saturation;
sum_pressure_watervolume += water_volume * water_pressure_reservoir;
sum_watervolume += water_volume;
cell_pressure[idx] = water_pressure_reservoir;
}
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&sum_pressure_watervolume, 1);
comm.sum(&sum_watervolume, 1);
// Ensure all processes have same notion of the aquifer cells' pressure values.
comm.sum(cell_pressure.data(), cell_pressure.size());
return sum_pressure_watervolume / sum_watervolume;
}

View File

@@ -83,14 +83,12 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::beginEpisode()
{
}
{}
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::beginIteration()
{
}
{}
template <typename TypeTag>
void
@@ -131,8 +129,7 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::endIteration()
{
}
{}
template <typename TypeTag>
void
@@ -155,11 +152,11 @@ BlackoilAquiferModel<TypeTag>::endTimeStep()
}
}
}
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::endEpisode()
{
}
{}
template <typename TypeTag>
template <class Restarter>
@@ -215,12 +212,14 @@ BlackoilAquiferModel<TypeTag>::init()
}
}
}
template <typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferCarterTracyActive() const
{
return !aquifers_CarterTracy.empty();
}
template <typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
@@ -236,7 +235,8 @@ BlackoilAquiferModel<TypeTag>::aquiferNumericalActive() const
}
template<typename TypeTag>
data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const
{
data::Aquifers data;
if (this->aquiferCarterTracyActive()) {
for (const auto& aqu : this->aquifers_CarterTracy) {

View File

@@ -226,6 +226,11 @@ HANDLE_AS_POD(data::NodeData)
HANDLE_AS_POD(data::Rates)
HANDLE_AS_POD(data::Segment)
std::size_t packSize(const data::NumericAquiferData& data, Dune::MPIHelper::MPICommunicator comm)
{
return packSize(data.initPressure, comm);
}
std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunicator comm)
{
const auto type = 0ull;
@@ -238,11 +243,20 @@ std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunic
+ packSize(data.datumDepth, comm)
+ packSize(type, comm);
if (data.aquFet != nullptr) {
return base + packSize(*data.aquFet, comm);
if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
aquFet != nullptr)
{
return base + packSize(*aquFet, comm);
}
else if (data.aquCT != nullptr) {
return base + packSize(*data.aquCT, comm);
else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
aquCT != nullptr)
{
return base + packSize(*aquCT, comm);
}
else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
aquNum != nullptr)
{
return base + packSize(*aquNum, comm);
}
return base;
@@ -513,12 +527,19 @@ void pack(const std::unordered_map<T1,T2,H,P,A>& data, std::vector<char>& buffer
}
}
void pack(const data::NumericAquiferData& data, std::vector<char>& buffer, int& position,
Dune::MPIHelper::MPICommunicator comm)
{
pack(data.initPressure, buffer, position, comm);
}
void pack(const data::AquiferData& data, std::vector<char>& buffer, int& position,
Dune::MPIHelper::MPICommunicator comm)
{
const auto type =
(data.aquFet != nullptr)*(1ull << 0)
+ (data.aquCT != nullptr)*(1ull << 1);
(data.typeData.is<data::AquiferType::Fetkovich>() * (1ull << 0))
+ (data.typeData.is<data::AquiferType::CarterTracy>() * (1ull << 1))
+ (data.typeData.is<data::AquiferType::Numerical>() * (1ull << 2));
pack(data.aquiferID, buffer, position, comm);
pack(data.pressure, buffer, position, comm);
@@ -528,11 +549,20 @@ void pack(const data::AquiferData& data, std::vector<char>& buffer, int& positio
pack(data.datumDepth, buffer, position, comm);
pack(type, buffer, position, comm);
if (data.aquFet != nullptr) {
pack(*data.aquFet, buffer, position, comm);
if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
aquFet != nullptr)
{
pack(*aquFet, buffer, position, comm);
}
else if (data.aquCT != nullptr) {
pack(*data.aquCT, buffer, position, comm);
else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
aquCT != nullptr)
{
pack(*aquCT, buffer, position, comm);
}
else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
aquNum != nullptr)
{
pack(*aquNum, buffer, position, comm);
}
}
@@ -865,6 +895,12 @@ void unpack(data::Well& data, std::vector<char>& buffer, int& position,
unpack(data.guide_rates, buffer, position, comm);
}
void unpack(data::NumericAquiferData& data, std::vector<char>& buffer, int& position,
Dune::MPIHelper::MPICommunicator comm)
{
unpack(data.initPressure, buffer, position, comm);
}
void unpack(data::AquiferData& data, std::vector<char>& buffer, int& position,
Dune::MPIHelper::MPICommunicator comm)
{
@@ -878,15 +914,17 @@ void unpack(data::AquiferData& data, std::vector<char>& buffer, int& position,
unpack(data.datumDepth, buffer, position, comm);
unpack(type, buffer, position, comm);
if (type == 1ull) {
data.type = data::AquiferType::Fetkovich;
data.aquFet = std::make_shared<data::FetkovichData>();
unpack(*data.aquFet, buffer, position, comm);
if (type == (1ull << 0)) {
auto* aquFet = data.typeData.create<data::AquiferType::Fetkovich>();
unpack(*aquFet, buffer, position, comm);
}
else if (type == 2ull) {
data.type = data::AquiferType::CarterTracy;
data.aquCT = std::make_shared<data::CarterTracyData>();
unpack(*data.aquCT, buffer, position, comm);
else if (type == (1ull << 1)) {
auto* aquCT = data.typeData.create<data::AquiferType::CarterTracy>();
unpack(*aquCT, buffer, position, comm);
}
else if (type == (1ull << 2)) {
auto* aquNum = data.typeData.create<data::AquiferType::Numerical>();
unpack(*aquNum, buffer, position, comm);
}
}

View File

@@ -61,6 +61,7 @@ struct GroupData;
struct GroupGuideRates;
class GuideRateValue;
struct NodeData;
struct NumericAquiferData;
class Rates;
struct Segment;
class Solution;
@@ -349,6 +350,7 @@ ADD_PACK_PROTOTYPES(data::GroupGuideRates)
ADD_PACK_PROTOTYPES(data::GroupData)
ADD_PACK_PROTOTYPES(data::NodeData)
ADD_PACK_PROTOTYPES(data::GroupAndNetworkValues)
ADD_PACK_PROTOTYPES(data::NumericAquiferData)
ADD_PACK_PROTOTYPES(data::Well)
ADD_PACK_PROTOTYPES(data::Wells)
ADD_PACK_PROTOTYPES(RestartKey)

View File

@@ -265,14 +265,14 @@ Opm::data::NodeData getNodeData()
Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1)
{
auto aquifer = Opm::data::AquiferData {
aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::Fetkovich
aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
};
aquifer.aquFet = std::make_shared<Opm::data::FetkovichData>();
auto* aquFet = aquifer.typeData.create<Opm::data::AquiferType::Fetkovich>();
aquifer.aquFet->initVolume = 1.23;
aquifer.aquFet->prodIndex = 45.67;
aquifer.aquFet->timeConstant = 890.123;
aquFet->initVolume = 1.23;
aquFet->prodIndex = 45.67;
aquFet->timeConstant = 890.123;
return aquifer;
}
@@ -280,17 +280,33 @@ Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1)
Opm::data::AquiferData getCarterTracyAquifer(const int aquiferID = 5)
{
auto aquifer = Opm::data::AquiferData {
aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::CarterTracy
aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
};
aquifer.aquCT = std::make_shared<Opm::data::CarterTracyData>();
auto* aquCT = aquifer.typeData.create<Opm::data::AquiferType::CarterTracy>();
aquifer.aquCT->timeConstant = 987.65;
aquifer.aquCT->influxConstant = 43.21;
aquifer.aquCT->waterDensity = 1014.5;
aquifer.aquCT->waterViscosity = 0.00318;
aquifer.aquCT->dimensionless_time = 42.0;
aquifer.aquCT->dimensionless_pressure = 2.34;
aquCT->timeConstant = 987.65;
aquCT->influxConstant = 43.21;
aquCT->waterDensity = 1014.5;
aquCT->waterViscosity = 0.00318;
aquCT->dimensionless_time = 42.0;
aquCT->dimensionless_pressure = 2.34;
return aquifer;
}
Opm::data::AquiferData getNumericalAquifer(const int aquiferID = 2)
{
auto aquifer = Opm::data::AquiferData {
aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
};
auto* aquNum = aquifer.typeData.create<Opm::data::AquiferType::Numerical>();
aquNum->initPressure.push_back(1.234);
aquNum->initPressure.push_back(2.345);
aquNum->initPressure.push_back(3.4);
aquNum->initPressure.push_back(9.876);
return aquifer;
}
@@ -363,12 +379,21 @@ BOOST_AUTO_TEST_CASE(dataCarterTracyData)
DO_CHECKS(data::CarterTracyData)
}
BOOST_AUTO_TEST_CASE(dataNumericAquiferData)
{
const auto val1 = getNumericalAquifer();
const auto val2 = PackUnpack(val1);
DO_CHECKS(data::NumericAquiferData)
}
BOOST_AUTO_TEST_CASE(dataAquifers)
{
const auto val1 = Opm::data::Aquifers {
{ 1, getFetkovichAquifer(1) },
{ 4, getFetkovichAquifer(4) },
{ 5, getCarterTracyAquifer(5) },
{ 9, getNumericalAquifer(9) },
};
const auto val2 = PackUnpack(val1);