Load Analytic Aquifer Objects from Restart File

This commit adds support for creating the analytic aquifer objects

    Opm::AquiferCT
    Opm::Aquifetp
    Opm::Aquancon
    Opm::AquiferConfig

from information stored in the restart vectors

    {I,S,X}AAQ
    {I,S}CAQ

We add a new helper class

    Opm::RestartIO::RstAquifer

which contain the same data members as the '*_data' structures of
the analytic aquifer objects.  Those analytic aquifer objects then
get friendship from the '*_data' structures in order to assign the
private members from the corresponding restart information.  We
finally add a gateway to EclipseState that consumes an RstAquifer
instance and overwrites the internal AquiferConfig object when the
restarted run contains analytic aquifers.

Update RstState constructor API to meet requirements of RstAquifer,
notably by adding an EclipseGrid parameter.  That in turn is needed
by the RstAquifer to translate connection (I,J,K) tuples to active
cell IDs.

Note that if an analytic aquifer does not have any connections then
this facility will not load said aquifer.  That may be a problem
when plotting summary curves, but we will address the issue later if
it comes up.
This commit is contained in:
Bård Skaflestad
2021-06-22 17:12:45 +02:00
parent d1b1cebdf5
commit 0d35d64aea
17 changed files with 1070 additions and 33 deletions

View File

@@ -284,6 +284,7 @@ if(ENABLE_ECL_OUTPUT)
src/opm/io/eclipse/ExtSmryOutput.cpp
src/opm/io/eclipse/RestartFileView.cpp
src/opm/io/eclipse/SummaryNode.cpp
src/opm/io/eclipse/rst/aquifer.cpp
src/opm/io/eclipse/rst/connection.cpp
src/opm/io/eclipse/rst/group.cpp
src/opm/io/eclipse/rst/header.cpp
@@ -864,6 +865,7 @@ if(ENABLE_ECL_OUTPUT)
opm/io/eclipse/ExtSmryOutput.hpp
opm/io/eclipse/RestartFileView.hpp
opm/io/eclipse/SummaryNode.hpp
opm/io/eclipse/rst/aquifer.hpp
opm/io/eclipse/rst/connection.hpp
opm/io/eclipse/rst/group.hpp
opm/io/eclipse/rst/header.hpp

View File

@@ -57,6 +57,8 @@ public:
getKeyword(const std::string& vector, const int occurrence = 0) const;
const std::vector<int>& intehead() const;
const std::vector<bool>& logihead() const;
const std::vector<double>& doubhead() const;
private:
class Implementation;

View File

@@ -0,0 +1,132 @@
/*
Copyright 2021 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 OPM_RESTART_AQUIFER_HPP
#define OPM_RESTART_AQUIFER_HPP
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <cstddef>
#include <memory>
#include <optional>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm {
class AquiferConfig;
class EclipseGrid;
class UnitSystem;
} // Opm
namespace Opm { namespace EclIO {
class RestartFileView;
}} // Opm::EclIO
namespace Opm { namespace RestartIO {
class RstAquifer
{
public:
struct CarterTracy {
int aquiferID{};
int inftableID{};
int pvttableID{};
double porosity{};
double datum_depth{};
double total_compr{};
double inner_radius{};
double permeability{};
double thickness{};
double angle_fraction{};
double initial_pressure{};
double time_constant{};
double influx_constant{};
double water_density{};
double water_viscosity{};
};
struct Fetkovich {
int aquiferID{};
int pvttableID{};
double prod_index{};
double total_compr{};
double initial_watvolume{};
double datum_depth{};
double initial_pressure{};
double time_constant{};
};
class Connections {
public:
struct Cell {
std::size_t global_index;
double influx_coeff;
double effective_facearea;
FaceDir::DirEnum face_dir;
};
const std::vector<Cell>& cells() const
{
return this->cells_;
}
void reserve(const std::vector<Cell>::size_type cpty)
{
this->cells_.reserve(cpty);
}
template <typename... Args>
void emplace_back(Args&&... args)
{
this->cells_.push_back(Cell { std::forward<Args>(args)... });
}
private:
std::vector<Cell> cells_{};
};
explicit RstAquifer(std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid* grid,
const UnitSystem& usys);
RstAquifer(const RstAquifer& rhs);
RstAquifer(RstAquifer&& rhs);
RstAquifer& operator=(const RstAquifer& rhs);
RstAquifer& operator=(RstAquifer&& rhs);
~RstAquifer();
bool hasAnalyticAquifers() const;
const std::vector<CarterTracy>& carterTracy() const;
const std::vector<Fetkovich>& fetkovich() const;
const std::unordered_map<int, Connections>& connections() const;
private:
class Implementation;
std::unique_ptr<Implementation> pImpl_;
};
}} // Opm::RestartIO
#endif // OPM_RESTART_AQUIFER_HPP

View File

@@ -24,6 +24,7 @@
#include <vector>
#include <opm/io/eclipse/rst/header.hpp>
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/io/eclipse/rst/group.hpp>
#include <opm/io/eclipse/rst/well.hpp>
#include <opm/io/eclipse/rst/udq.hpp>
@@ -32,6 +33,9 @@
#include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
namespace Opm {
class EclipseGrid;
} // namespace Opm
namespace Opm { namespace EclIO {
class RestartFileView;
@@ -40,17 +44,17 @@ namespace Opm { namespace EclIO {
namespace Opm { namespace RestartIO {
struct RstState {
RstState(const ::Opm::UnitSystem& unit_system,
const std::vector<int>& intehead,
const std::vector<bool>& logihead,
const std::vector<double>& doubhead);
RstState(std::shared_ptr<EclIO::RestartFileView> rstView,
const ::Opm::EclipseGrid* grid);
static RstState load(std::shared_ptr<EclIO::RestartFileView> rstView);
static RstState load(std::shared_ptr<EclIO::RestartFileView> rstView,
const ::Opm::EclipseGrid* grid = nullptr);
const RstWell& get_well(const std::string& wname) const;
const ::Opm::UnitSystem unit_system;
::Opm::UnitSystem unit_system;
RstHeader header;
RstAquifer aquifers;
std::vector<RstWell> wells;
std::vector<RstGroup> groups;
std::vector<RstUDQ> udqs;

View File

@@ -37,6 +37,10 @@ namespace Opm {
class Deck;
}
namespace Opm { namespace RestartIO {
class RstAquifer;
}} // Opm::RestartIO
namespace Opm {
class Aquancon {
@@ -82,11 +86,12 @@ namespace Opm {
}
};
Aquancon() = default;
Aquancon(const EclipseGrid& grid, const Deck& deck);
explicit Aquancon(const std::unordered_map<int, std::vector<Aquancon::AquancCell>>& data);
void loadFromRestart(const RestartIO::RstAquifer& rst_aquifers);
static Aquancon serializeObject();
const std::unordered_map<int, std::vector<Aquancon::AquancCell>>& data() const;

View File

@@ -36,6 +36,10 @@ namespace Opm {
class TableManager;
}
namespace Opm { namespace RestartIO {
class RstAquifer;
}} // Opm::RestartIO
namespace Opm {
class AquiferCT {
@@ -43,6 +47,8 @@ namespace Opm {
struct AQUCT_data
{
friend class AquiferCT;
AQUCT_data() = default;
AQUCT_data(const DeckRecord& record, const TableManager& tables);
AQUCT_data(const int aqID,
@@ -117,6 +123,9 @@ namespace Opm {
AquiferCT(const TableManager& tables, const Deck& deck);
AquiferCT(const std::vector<AquiferCT::AQUCT_data>& data);
void loadFromRestart(const RestartIO::RstAquifer& rst,
const TableManager& tables);
static AquiferCT serializeObject();
std::size_t size() const;

View File

@@ -26,11 +26,17 @@
#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp>
namespace Opm {
class TableManager;
class EclipseGrid;
class Deck;
class FieldPropsManager;
} // namespace Opm
class TableManager;
class EclipseGrid;
class Deck;
class FieldPropsManager;
namespace Opm { namespace RestartIO {
class RstAquifer;
}} // namespace Opm::RestartIO
namespace Opm {
class AquiferConfig {
public:
@@ -40,6 +46,8 @@ public:
AquiferConfig(const Aquifetp& fetp, const AquiferCT& ct, const Aquancon& conn);
void load_connections(const Deck& deck, const EclipseGrid& grid);
void loadFromRestart(const RestartIO::RstAquifer& aquifers,
const TableManager& tables);
static AquiferConfig serializeObject();

View File

@@ -36,6 +36,10 @@ namespace Opm {
class TableManager;
}
namespace Opm { namespace RestartIO {
class RstAquifer;
}} // Opm::RestartIO
namespace Opm {
class Aquifetp {
@@ -43,6 +47,8 @@ class Aquifetp {
struct AQUFETP_data
{
friend class Aquifetp;
AQUFETP_data() = default;
AQUFETP_data(const DeckRecord& record, const TableManager& tables);
AQUFETP_data(const int aquiferID_,
@@ -98,6 +104,9 @@ class Aquifetp {
Aquifetp(const TableManager& tables, const Deck& deck);
explicit Aquifetp(const std::vector<Aquifetp::AQUFETP_data>& data);
void loadFromRestart(const RestartIO::RstAquifer& rst,
const TableManager& tables);
static Aquifetp serializeObject();
const std::vector<Aquifetp::AQUFETP_data>& data() const;

View File

@@ -15,7 +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/>.
*/
*/
#ifndef OPM_ECLIPSE_STATE_HPP
#define OPM_ECLIPSE_STATE_HPP
@@ -37,13 +37,18 @@
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
namespace Opm {
class Deck;
class DeckKeyword;
class InitConfig;
class IOConfig;
class DeckSection;
} // namespace Opm
namespace Opm { namespace RestartIO {
class RstAquifer;
}} // namespace Opm::RestartIO
namespace Opm {
class EclipseState {
public:
@@ -100,6 +105,8 @@ namespace Opm {
const AquiferConfig& aquifer() const;
const TracerConfig& tracer() const;
void loadRestartAquifers(const RestartIO::RstAquifer& aquifers);
template<class Serializer>
void serializeOp(Serializer& serializer)
{
@@ -148,6 +155,6 @@ namespace Opm {
std::string m_title;
TracerConfig tracer_config;
};
}
} // namespace Opm
#endif // OPM_ECLIPSE_STATE_HPP

View File

@@ -132,6 +132,32 @@ public:
return this->getKeyword<int>(ihkw, 0);
}
const std::vector<bool>& logihead()
{
const auto lhkw = std::string { "LOGIHEAD" };
if (! this->hasKeyword<bool>(lhkw)) {
throw std::domain_error {
"Purported Restart File Does not Have Logical Header"
};
}
return this->getKeyword<bool>(lhkw, 0);
}
const std::vector<double>& doubhead()
{
const auto dhkw = std::string { "DOUBHEAD" };
if (! this->hasKeyword<double>(dhkw)) {
throw std::domain_error {
"Purported Restart File Does not Have Double Header"
};
}
return this->getKeyword<double>(dhkw, 0);
}
private:
using RstFile = std::shared_ptr<ERst>;
@@ -170,7 +196,6 @@ Implementation(std::shared_ptr<ERst> restart_file,
const auto& type = std::get<1>(vector);
switch (type) {
case ::Opm::EclIO::eclArrType::LOGI:
case ::Opm::EclIO::eclArrType::MESS:
// Currently ignored
continue;
@@ -240,6 +265,16 @@ const std::vector<int>& Opm::EclIO::RestartFileView::intehead() const
return this->pImpl_->intehead();
}
const std::vector<bool>& Opm::EclIO::RestartFileView::logihead() const
{
return this->pImpl_->logihead();
}
const std::vector<double>& Opm::EclIO::RestartFileView::doubhead() const
{
return this->pImpl_->doubhead();
}
template <typename ElmType>
bool Opm::EclIO::RestartFileView::hasKeyword(const std::string& vector) const
{

View File

@@ -0,0 +1,683 @@
/*
Copyright 2021 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/>.
*/
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/io/eclipse/RestartFileView.hpp>
#include <opm/output/eclipse/VectorItems/aquifer.hpp>
#include <opm/output/eclipse/VectorItems/intehead.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <cstddef>
#include <memory>
#include <stdexcept>
#include <unordered_map>
#include <utility>
#include <vector>
#include <fmt/format.h>
#include <boost/range.hpp>
namespace VI = ::Opm::RestartIO::Helpers::VectorItems;
namespace {
template <typename T>
boost::iterator_range<typename std::vector<T>::const_iterator>
getDataWindow(const std::vector<T>& arr,
const std::size_t windowSize,
const std::size_t entity,
const std::size_t subEntity = 0,
const std::size_t maxSubEntitiesPerEntity = 1)
{
const auto off =
windowSize * (subEntity + maxSubEntitiesPerEntity*entity);
auto begin = arr.begin() + off;
auto end = begin + windowSize;
return { begin, end };
}
// ---------------------------------------------------------------------
class ConnectionVectors
{
public:
template <typename T>
using Window = boost::iterator_range<
typename std::vector<T>::const_iterator
>;
explicit ConnectionVectors(const std::vector<int>& intehead,
std::shared_ptr<Opm::EclIO::RestartFileView> rst_view);
Window<int> icaq(const int occurrence,
const std::size_t connectionID) const;
Window<float> scaq(const int occurrence,
const std::size_t connectionID) const;
private:
std::size_t numIConnElem_;
std::size_t numSConnElem_;
std::shared_ptr<Opm::EclIO::RestartFileView> rstView_;
};
ConnectionVectors::ConnectionVectors(const std::vector<int>& intehead,
std::shared_ptr<Opm::EclIO::RestartFileView> rst_view)
: numIConnElem_(intehead[VI::intehead::NICAQZ])
, numSConnElem_(intehead[VI::intehead::NSCAQZ])
, rstView_ (std::move(rst_view))
{}
ConnectionVectors::Window<int>
ConnectionVectors::icaq(const int occurrence,
const std::size_t connectionID) const
{
return getDataWindow(this->rstView_->getKeyword<int>("ICAQ", occurrence),
this->numIConnElem_, connectionID);
}
ConnectionVectors::Window<float>
ConnectionVectors::scaq(const int occurrence,
const std::size_t connectionID) const
{
return getDataWindow(this->rstView_->getKeyword<float>("SCAQ", occurrence),
this->numSConnElem_, connectionID);
}
// ---------------------------------------------------------------------
class AquiferVectors
{
public:
template <typename T>
using Window = boost::iterator_range<
typename std::vector<T>::const_iterator
>;
explicit AquiferVectors(const std::vector<int>& intehead,
std::shared_ptr<Opm::EclIO::RestartFileView> rst_view);
int maxAquiferID() const;
Window<int> iaaq(const std::size_t aquiferID) const;
Window<float> saaq(const std::size_t aquiferID) const;
Window<double> xaaq(const std::size_t aquiferID) const;
private:
int maxAquiferID_{};
std::size_t numIAquifElem_;
std::size_t numSAquifElem_;
std::size_t numXAquifElem_;
std::shared_ptr<Opm::EclIO::RestartFileView> rstView_;
};
AquiferVectors::AquiferVectors(const std::vector<int>& intehead,
std::shared_ptr<Opm::EclIO::RestartFileView> rst_view)
: maxAquiferID_ (intehead[VI::intehead::MAX_AN_AQUIFER_ID])
, numIAquifElem_(intehead[VI::intehead::NIAAQZ])
, numSAquifElem_(intehead[VI::intehead::NSAAQZ])
, numXAquifElem_(intehead[VI::intehead::NXAAQZ])
, rstView_ (std::move(rst_view))
{}
int AquiferVectors::maxAquiferID() const
{
return this->maxAquiferID_;
}
AquiferVectors::Window<int>
AquiferVectors::iaaq(const std::size_t aquiferID) const
{
return getDataWindow(this->rstView_->getKeyword<int>("IAAQ"),
this->numIAquifElem_, aquiferID);
}
AquiferVectors::Window<float>
AquiferVectors::saaq(const std::size_t aquiferID) const
{
return getDataWindow(this->rstView_->getKeyword<float>("SAAQ"),
this->numSAquifElem_, aquiferID);
}
AquiferVectors::Window<double>
AquiferVectors::xaaq(const std::size_t aquiferID) const
{
return getDataWindow(this->rstView_->getKeyword<double>("XAAQ"),
this->numXAquifElem_, aquiferID);
}
// ---------------------------------------------------------------------
class ConnectionOccurrence
{
public:
explicit ConnectionOccurrence(const int numAquifers,
Opm::EclIO::RestartFileView& rst_view);
int icaq(const int aquiferID) const
{
return this->occurrence_[aquiferID].I;
}
int scaq(const int aquiferID) const
{
return this->occurrence_[aquiferID].S;
}
private:
struct Occurrence
{
int I{-1};
int S{-1};
};
std::vector<Occurrence> occurrence_{};
};
ConnectionOccurrence::ConnectionOccurrence(const int maxAquiferID,
Opm::EclIO::RestartFileView& rst_view)
: occurrence_(maxAquiferID)
{
const auto integerAquifID = std::string{ "ICAQNUM" };
const auto floatAquifID = std::string{ "SCAQNUM" };
if ((rst_view.occurrenceCount(integerAquifID) < maxAquiferID) ||
(rst_view.occurrenceCount(floatAquifID) < maxAquiferID))
{
throw std::invalid_argument {
"Structural inconsistency for analytic "
"aquifer connections in restart file"
};
}
// Note:
// One separate xCAQNUM table for each aquifer ID in 1..maxAquiferID.
//
// On the one hand, Flow's output system at the time of writing
// guarantees that
//
// occurrence_[i].I == occurrence_[i].S == i
//
// for all 'i' in 0..maxAquiferID-1. On the other hand, the file
// format allows for general ordering. This mapping facility exists
// mainly to handle the general case.
for (auto occurrence = 0*maxAquiferID; occurrence < maxAquiferID; ++occurrence) {
const auto& iAquifID = rst_view.getKeyword<int>(integerAquifID, occurrence);
const auto& sAquifID = rst_view.getKeyword<int>(floatAquifID, occurrence);
this->occurrence_[iAquifID.front() - 1].I = occurrence;
this->occurrence_[sAquifID.front() - 1].S = occurrence;
}
}
// ---------------------------------------------------------------------
Opm::FaceDir::DirEnum face_direction(const int directionValue)
{
using FDValue = VI::IAnalyticAquiferConn::Value::FaceDirection;
using FDir = Opm::FaceDir::DirEnum;
switch (directionValue) {
case FDValue::IMinus: return FDir::XMinus;
case FDValue::IPlus: return FDir::XPlus;
case FDValue::JMinus: return FDir::YMinus;
case FDValue::JPlus: return FDir::YPlus;
case FDValue::KMinus: return FDir::ZMinus;
case FDValue::KPlus: return FDir::ZPlus;
}
throw std::invalid_argument {
fmt::format("Unknown Face Direction {}", directionValue)
};
}
template <typename IcaqArray>
Opm::FaceDir::DirEnum face_direction(const IcaqArray& icaq)
{
using Ix = VI::IAnalyticAquiferConn::index;
return face_direction(icaq[Ix::FaceDirection]);
}
template <typename IcaqArray>
std::size_t identify_global_cell(const Opm::EclipseGrid& grid,
const IcaqArray& icaq)
{
using Ix = VI::IAnalyticAquiferConn::index;
const auto i = static_cast<std::size_t>(icaq[Ix::Index_I] - 1);
const auto j = static_cast<std::size_t>(icaq[Ix::Index_J] - 1);
const auto k = static_cast<std::size_t>(icaq[Ix::Index_K] - 1);
return grid.getGlobalIndex(i, j, k);
}
template <typename ScaqArray>
double influx_coefficient(const ScaqArray scaq, const double tot_influx)
{
using Ix = VI::SAnalyticAquiferConn::index;
return scaq[Ix::InfluxFraction] * tot_influx;
}
template <typename ScaqArray>
double effective_facearea(const ScaqArray scaq, const double tot_influx)
{
using Ix = VI::SAnalyticAquiferConn::index;
return scaq[Ix::FaceAreaToInfluxCoeff] * tot_influx;
}
template <typename IcaqArray, typename ScaqArray>
void load_analytic_aquifer_cell(const IcaqArray icaq,
const ScaqArray scaq,
const Opm::EclipseGrid& grid,
const double tot_influx,
Opm::RestartIO::RstAquifer::Connections& connections)
{
const auto glob_cell = identify_global_cell(grid, icaq);
const auto influx_coeff = influx_coefficient(scaq, tot_influx);
const auto eff_facearea = effective_facearea(scaq, tot_influx);
const auto face_dir = face_direction(icaq);
connections.emplace_back(glob_cell, influx_coeff, eff_facearea, face_dir);
}
Opm::RestartIO::RstAquifer::Connections
load_analytic_aquifer_cells(const ConnectionOccurrence& occurence,
const ConnectionVectors& connections,
const Opm::EclipseGrid& grid,
const std::size_t num_conn,
const double tot_influx,
const int aquifer_id)
{
auto cells = Opm::RestartIO::RstAquifer::Connections{};
cells.reserve(num_conn);
const auto occur_icaq = occurence.icaq(aquifer_id);
const auto occur_scaq = occurence.scaq(aquifer_id);
for (auto connID = 0*num_conn; connID < num_conn; ++connID) {
auto icaq = connections.icaq(occur_icaq, connID);
auto scaq = connections.scaq(occur_scaq, connID);
load_analytic_aquifer_cell(icaq, scaq, grid, tot_influx, cells);
}
return cells;
}
// ---------------------------------------------------------------------
Opm::RestartIO::RstAquifer::CarterTracy
load_carter_tracy(const int aquiferID,
const AquiferVectors& aquifers,
const Opm::UnitSystem& usys)
{
auto aquifer = Opm::RestartIO::RstAquifer::CarterTracy{};
using M = Opm::UnitSystem::measure;
using IntIX = VI::IAnalyticAquifer::index;
using RealIX = VI::SAnalyticAquifer::index;
using DoubIX = VI::XAnalyticAquifer::index;
const auto iaaq = aquifers.iaaq(aquiferID);
const auto saaq = aquifers.saaq(aquiferID);
const auto xaaq = aquifers.xaaq(aquiferID);
aquifer.aquiferID = aquiferID + 1;
aquifer.inftableID = iaaq[IntIX::CTInfluenceFunction];
aquifer.pvttableID = iaaq[IntIX::WatPropTable];
aquifer.porosity = usys.to_si(M::identity, saaq[RealIX::CTPorosity]);
aquifer.datum_depth = usys.to_si(M::length, saaq[RealIX::DatumDepth]);
// Note: *from_si()* to work around the fact that we don't have a
// compressibility unit.
aquifer.total_compr =
usys.from_si(M::pressure, saaq[RealIX::Compressibility]);
aquifer.permeability =
usys.to_si(M::permeability, saaq[RealIX::CTPermeability]);
aquifer.inner_radius = usys.to_si(M::length, saaq[RealIX::CTRadius]);
aquifer.thickness = usys.to_si(M::length, saaq[RealIX::CTThickness]);
aquifer.angle_fraction = usys.to_si(M::identity, saaq[RealIX::CTAngle]);
aquifer.initial_pressure = usys.to_si(M::pressure, saaq[RealIX::InitPressure]);
aquifer.time_constant = usys.to_si(M::time, 1.0 / xaaq[DoubIX::CTRecipTimeConst]);
// Note: *from_si()* for the pressure unit here since 'beta' is total
// influx (volume) per unit pressure drop.
aquifer.influx_constant =
usys.to_si(M::volume, usys.from_si(M::pressure, xaaq[DoubIX::CTInfluxConstant]));
aquifer.water_density = usys.to_si(M::density, saaq[RealIX::CTWatMassDensity]);
aquifer.water_viscosity = usys.to_si(M::viscosity, saaq[RealIX::CTWatViscosity]);
return aquifer;
}
// ---------------------------------------------------------------------
Opm::RestartIO::RstAquifer::Fetkovich
load_fetkovich(const int aquiferID,
const AquiferVectors& aquifers,
const Opm::UnitSystem& usys)
{
auto aquifer = Opm::RestartIO::RstAquifer::Fetkovich{};
using M = Opm::UnitSystem::measure;
using IntIX = VI::IAnalyticAquifer::index;
using RealIX = VI::SAnalyticAquifer::index;
const auto iaaq = aquifers.iaaq(aquiferID);
const auto saaq = aquifers.saaq(aquiferID);
aquifer.aquiferID = aquiferID + 1;
aquifer.pvttableID = iaaq[IntIX::WatPropTable];
aquifer.prod_index =
usys.to_si(M::liquid_productivity_index, saaq[RealIX::FetProdIndex]);
// Note: *from_si()* to work around the fact that we don't have a
// compressibility unit.
aquifer.total_compr =
usys.from_si(M::pressure, saaq[RealIX::Compressibility]);
aquifer.initial_watvolume =
usys.to_si(M::liquid_surface_volume, saaq[RealIX::FetInitVol]);
aquifer.datum_depth = usys.to_si(M::length, saaq[RealIX::DatumDepth]);
aquifer.initial_pressure = usys.to_si(M::pressure, saaq[RealIX::InitPressure]);
aquifer.time_constant = usys.to_si(M::time, saaq[RealIX::FetTimeConstant]);
return aquifer;
}
// ---------------------------------------------------------------------
int num_aquifers(const std::vector<int>& intehead)
{
return intehead[VI::intehead::NAQUIF];
}
int num_aquifer_connections(const AquiferVectors& aquifers,
const std::size_t aquiferID)
{
using Ix = VI::IAnalyticAquifer::index;
auto iaaq = aquifers.iaaq(aquiferID);
return iaaq[Ix::NumAquiferConn];
}
std::unordered_map<int, Opm::RestartIO::RstAquifer::Connections>
load_aquifer_connections(const ConnectionOccurrence& occurence,
const AquiferVectors& aquifers,
const ConnectionVectors& connections,
const Opm::EclipseGrid& grid,
const Opm::UnitSystem& usys,
const int max_aquifer_id)
{
auto aqConn = std::unordered_map<int, Opm::RestartIO::RstAquifer::Connections>{};
auto tot_influx = [&aquifers, &usys](const std::size_t aquiferID) -> double
{
using M = Opm::UnitSystem::measure;
using Ix = VI::XAnalyticAquifer::index;
auto xaaq = aquifers.xaaq(aquiferID);
return usys.to_si(M::length, usys.to_si(M::length, xaaq[Ix::TotalInfluxCoeff]));
};
auto load = [&occurence, &connections, &grid]
(const std::size_t num_connections,
const double total_influx,
const int aquifer_id)
{
return load_analytic_aquifer_cells(occurence, connections, grid,
num_connections, total_influx, aquifer_id);
};
for (auto aquiferID = 0*max_aquifer_id; aquiferID < max_aquifer_id; ++aquiferID) {
const auto num_connections = num_aquifer_connections(aquifers, aquiferID);
if (num_connections == 0) {
continue;
}
const auto total_influx = tot_influx(aquiferID);
aqConn.insert_or_assign(aquiferID + 1, load(num_connections, total_influx, aquiferID));
}
return aqConn;
}
} // Anonymous
class Opm::RestartIO::RstAquifer::Implementation
{
public:
explicit Implementation(std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid* grid,
const UnitSystem& usys);
bool hasAnalyticAquifers() const
{
return ! (this->connections_.empty() &&
this->carterTracy_.empty() &&
this->fetkovich_ .empty());
}
const std::vector<RstAquifer::CarterTracy>& carterTracy() const
{
return this->carterTracy_;
}
const std::vector<RstAquifer::Fetkovich>& fetkovich() const
{
return this->fetkovich_;
}
const std::unordered_map<int, RstAquifer::Connections>& connections() const
{
return this->connections_;
}
private:
std::unordered_map<int, RstAquifer::Connections> connections_{};
std::vector<RstAquifer::CarterTracy> carterTracy_{};
std::vector<RstAquifer::Fetkovich> fetkovich_{};
void loadAnalyticAquiferConnections(const AquiferVectors& vectors,
std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid& grid,
const UnitSystem& usys);
void loadAnalyticAquifers(const AquiferVectors& aquifers,
const UnitSystem& usys);
void loadAnalyticAquifer(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys);
void loadCarterTracy(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys);
void loadFetkovich(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys);
};
Opm::RestartIO::RstAquifer::Implementation::
Implementation(std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid* grid,
const UnitSystem& usys)
{
const auto numAquifers = num_aquifers(rstView->intehead());
if ((numAquifers == 0) || (grid == nullptr)) {
return;
}
const auto aquifers = AquiferVectors { rstView->intehead(), rstView };
this->loadAnalyticAquiferConnections(aquifers, std::move(rstView), *grid, usys);
this->loadAnalyticAquifers(aquifers, usys);
}
void
Opm::RestartIO::RstAquifer::Implementation::
loadAnalyticAquiferConnections(const AquiferVectors& aquifers,
std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid& grid,
const UnitSystem& usys)
{
const auto occurence = ConnectionOccurrence { aquifers.maxAquiferID(), *rstView };
const auto connections = ConnectionVectors { rstView->intehead(), std::move(rstView) };
this->connections_ =
load_aquifer_connections(occurence, aquifers, connections,
grid, usys, aquifers.maxAquiferID());
}
void
Opm::RestartIO::RstAquifer::Implementation::
loadAnalyticAquifers(const AquiferVectors& aquifers,
const UnitSystem& usys)
{
const auto maxAquiferID = aquifers.maxAquiferID();
for (auto aquiferID = 0*maxAquiferID; aquiferID < maxAquiferID; ++aquiferID) {
if (num_aquifer_connections(aquifers, aquiferID) == 0) {
// Skip aquifers without connections. Likely to be sparsely
// allocated aquifer IDs (e.g., 1, 2, 5).
continue;
}
this->loadAnalyticAquifer(aquiferID, aquifers, usys);
}
}
void
Opm::RestartIO::RstAquifer::Implementation::
loadAnalyticAquifer(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys)
{
const auto iaaq = aquifers.iaaq(aquiferID);
const auto type = iaaq[VI::IAnalyticAquifer::index::TypeRelated1];
switch (type) {
case VI::IAnalyticAquifer::Value::ModelType::CarterTracy:
this->loadCarterTracy(aquiferID, aquifers, usys);
return;
case VI::IAnalyticAquifer::Value::ModelType::Fetkovich:
this->loadFetkovich(aquiferID, aquifers, usys);
return;
}
throw std::invalid_argument {
fmt::format("Analytic aquifer {} (type {}) is neither Carter-Tracy "
"nor Fetkovich in restart input file", aquiferID + 1, type)
};
}
void
Opm::RestartIO::RstAquifer::Implementation::
loadCarterTracy(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys)
{
this->carterTracy_.push_back(load_carter_tracy(aquiferID, aquifers, usys));
}
void
Opm::RestartIO::RstAquifer::Implementation::
loadFetkovich(const int aquiferID,
const AquiferVectors& aquifers,
const UnitSystem& usys)
{
this->fetkovich_.push_back(load_fetkovich(aquiferID, aquifers, usys));
}
// ---------------------------------------------------------------------
Opm::RestartIO::RstAquifer::RstAquifer(std::shared_ptr<EclIO::RestartFileView> rstView,
const EclipseGrid* grid,
const UnitSystem& usys)
: pImpl_{ new Implementation{ std::move(rstView), grid, usys } }
{}
Opm::RestartIO::RstAquifer::RstAquifer(const RstAquifer& rhs)
: pImpl_{ new Implementation{ *rhs.pImpl_ } }
{}
Opm::RestartIO::RstAquifer::RstAquifer(RstAquifer&& rhs)
: pImpl_{ std::move(rhs.pImpl_) }
{}
Opm::RestartIO::RstAquifer&
Opm::RestartIO::RstAquifer::operator=(const RstAquifer& rhs)
{
this->pImpl_.reset(new Implementation{ *rhs.pImpl_ });
return *this;
}
Opm::RestartIO::RstAquifer&
Opm::RestartIO::RstAquifer::operator=(RstAquifer&& rhs)
{
this->pImpl_ = std::move(rhs.pImpl_);
return *this;
}
Opm::RestartIO::RstAquifer::~RstAquifer()
{}
bool Opm::RestartIO::RstAquifer::hasAnalyticAquifers() const
{
return this->pImpl_->hasAnalyticAquifers();
}
const std::vector<Opm::RestartIO::RstAquifer::CarterTracy>&
Opm::RestartIO::RstAquifer::carterTracy() const
{
return this->pImpl_->carterTracy();
}
const std::vector<Opm::RestartIO::RstAquifer::Fetkovich>&
Opm::RestartIO::RstAquifer::fetkovich() const
{
return this->pImpl_->fetkovich();
}
const std::unordered_map<int, Opm::RestartIO::RstAquifer::Connections>&
Opm::RestartIO::RstAquifer::connections() const
{
return this->pImpl_->connections();
}

View File

@@ -20,6 +20,7 @@
#include <iterator>
#include <memory>
#include <numeric>
#include <optional>
#include <opm/io/eclipse/RestartFileView.hpp>
@@ -64,14 +65,13 @@ namespace VI = ::Opm::RestartIO::Helpers::VectorItems;
namespace Opm { namespace RestartIO {
RstState::RstState(const ::Opm::UnitSystem& unit_system_,
const std::vector<int>& intehead,
const std::vector<bool>& logihead,
const std::vector<double>& doubhead):
unit_system(unit_system_),
header(unit_system_, intehead, logihead, doubhead)
RstState::RstState(std::shared_ptr<EclIO::RestartFileView> rstView,
const ::Opm::EclipseGrid* grid)
: unit_system(rstView->intehead()[VI::intehead::UNIT])
, header(unit_system, rstView->intehead(), rstView->logihead(), rstView->doubhead())
, aquifers(rstView, grid, unit_system)
{
this->load_tuning(intehead, doubhead);
this->load_tuning(rstView->intehead(), rstView->doubhead());
}
@@ -273,15 +273,11 @@ const RstWell& RstState::get_well(const std::string& wname) const {
return *well_iter;
}
RstState RstState::load(std::shared_ptr<EclIO::RestartFileView> rstView) {
const auto& intehead = rstView->getKeyword<int>("INTEHEAD");
const auto& logihead = rstView->getKeyword<bool>("LOGIHEAD");
const auto& doubhead = rstView->getKeyword<double>("DOUBHEAD");
RstState RstState::load(std::shared_ptr<EclIO::RestartFileView> rstView,
const ::Opm::EclipseGrid* grid)
{
RstState state(rstView, grid);
auto unit_id = intehead[VI::intehead::UNIT];
::Opm::UnitSystem unit_system(unit_id);
RstState state(unit_system, intehead, logihead, doubhead);
if (state.header.ngroup > 0) {
const auto& zgrp = rstView->getKeyword<std::string>("ZGRP");
const auto& igrp = rstView->getKeyword<int>("IGRP");

View File

@@ -15,10 +15,12 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
*/
#include <opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp>
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
@@ -46,6 +48,36 @@
#include "AquiferHelpers.hpp"
namespace {
Opm::Aquancon::AquancCell
makeAquiferCell(const int aquiferID,
const Opm::RestartIO::RstAquifer::Connections::Cell& rst_cell)
{
return {
aquiferID,
rst_cell.global_index,
rst_cell.influx_coeff,
rst_cell.effective_facearea,
rst_cell.face_dir
};
}
std::vector<Opm::Aquancon::AquancCell>
makeAquiferConnections(const int aquiferID,
const Opm::RestartIO::RstAquifer::Connections& rst_connections)
{
auto connections = std::vector<Opm::Aquancon::AquancCell>{};
const auto& rst_cells = rst_connections.cells();
connections.reserve(rst_cells.size());
for (const auto& rst_cell : rst_cells) {
connections.push_back(makeAquiferCell(aquiferID, rst_cell));
}
return connections;
}
} // Anonymous
namespace Opm {
namespace{
@@ -183,6 +215,15 @@ namespace Opm {
cells(data)
{}
void Aquancon::loadFromRestart(const RestartIO::RstAquifer& rst_aquifers)
{
this->cells.clear();
for (const auto& [aquiferID, rst_connections] : rst_aquifers.connections()) {
this->cells.insert_or_assign(aquiferID, makeAquiferConnections(aquiferID, rst_connections));
}
}
const std::unordered_map<int, std::vector<Aquancon::AquancCell>>& Aquancon::data() const {
return this->cells;
}

View File

@@ -15,10 +15,12 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
*/
#include <opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp>
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Aqudims.hpp>
@@ -42,6 +44,29 @@
#include <utility>
#include <vector>
namespace {
Opm::AquiferCT::AQUCT_data
makeAquifer(const Opm::RestartIO::RstAquifer::CarterTracy& rst_aquifer)
{
auto aquifer = Opm::AquiferCT::AQUCT_data{};
aquifer.aquiferID = rst_aquifer.aquiferID;
aquifer.inftableID = rst_aquifer.inftableID;
aquifer.pvttableID = rst_aquifer.pvttableID;
aquifer.porosity = rst_aquifer.porosity;
aquifer.datum_depth = rst_aquifer.datum_depth;
aquifer.total_compr = rst_aquifer.total_compr;
aquifer.inner_radius = rst_aquifer.inner_radius;
aquifer.permeability = rst_aquifer.permeability;
aquifer.thickness = rst_aquifer.thickness;
aquifer.angle_fraction = rst_aquifer.angle_fraction;
aquifer.initial_pressure = rst_aquifer.initial_pressure;
return aquifer;
}
}
namespace Opm {
namespace {
@@ -199,6 +224,30 @@ AquiferCT::AquiferCT(const std::vector<AquiferCT::AQUCT_data>& data) :
m_aquct(data)
{}
void AquiferCT::loadFromRestart(const RestartIO::RstAquifer& rst,
const TableManager& tables)
{
this->m_aquct.clear();
const auto& rst_aquifers = rst.carterTracy();
if (! rst_aquifers.empty()) {
this->m_aquct.reserve(rst_aquifers.size());
}
for (const auto& rst_aquifer : rst_aquifers) {
this->m_aquct.push_back(makeAquifer(rst_aquifer));
auto& new_aquifer = this->m_aquct.back();
new_aquifer.finishInitialisation(tables);
// Assign 'private:' members in AQUCT_data.
new_aquifer.time_constant_ = rst_aquifer.time_constant;
new_aquifer.influx_constant_ = rst_aquifer.influx_constant;
new_aquifer.water_density_ = rst_aquifer.water_density;
new_aquifer.water_viscosity_ = rst_aquifer.water_viscosity;
}
}
AquiferCT AquiferCT::serializeObject()
{
return { { AQUCT_data::serializeObject() } };

View File

@@ -47,6 +47,14 @@ void AquiferConfig::load_connections(const Deck& deck, const EclipseGrid& grid)
this->aqconn = Aquancon(grid, deck);
}
void AquiferConfig::loadFromRestart(const RestartIO::RstAquifer& aquifers,
const TableManager& tables)
{
this->aquifetp.loadFromRestart(aquifers, tables);
this->aquiferct.loadFromRestart(aquifers, tables);
this->aqconn.loadFromRestart(aquifers);
}
AquiferConfig AquiferConfig::serializeObject()
{
AquiferConfig result;

View File

@@ -19,6 +19,8 @@
#include <opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp>
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/FlatTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
@@ -32,6 +34,25 @@
#include <vector>
namespace {
Opm::Aquifetp::AQUFETP_data
makeAquifer(const Opm::RestartIO::RstAquifer::Fetkovich& rst_aquifer)
{
auto aquifer = Opm::Aquifetp::AQUFETP_data{};
aquifer.aquiferID = rst_aquifer.aquiferID;
aquifer.pvttableID = rst_aquifer.pvttableID;
aquifer.prod_index = rst_aquifer.prod_index;
aquifer.total_compr = rst_aquifer.total_compr;
aquifer.initial_watvolume = rst_aquifer.initial_watvolume;
aquifer.datum_depth = rst_aquifer.datum_depth;
aquifer.initial_pressure = rst_aquifer.initial_pressure;
return aquifer;
}
}
namespace Opm {
using AQUFETP = ParserKeywords::AQUFETP;
@@ -131,6 +152,26 @@ Aquifetp::Aquifetp(const std::vector<Aquifetp::AQUFETP_data>& data) :
m_aqufetp(data)
{}
void Aquifetp::loadFromRestart(const RestartIO::RstAquifer& rst,
const TableManager& tables)
{
this->m_aqufetp.clear();
const auto& rst_aquifers = rst.fetkovich();
if (! rst_aquifers.empty()) {
this->m_aqufetp.reserve(rst_aquifers.size());
}
for (const auto& rst_aquifer : rst_aquifers) {
this->m_aqufetp.push_back(makeAquifer(rst_aquifer));
auto& new_aquifer = this->m_aqufetp.back();
new_aquifer.finishInitialisation(tables);
// Assign 'private:' members in AQUFETP_data.
new_aquifer.time_constant_ = rst_aquifer.time_constant;
}
}
Aquifetp Aquifetp::serializeObject()
{

View File

@@ -26,6 +26,8 @@
#include <opm/common/OpmLog/LogUtil.hpp>
#include <opm/common/utility/OpmInputError.hpp>
#include <opm/io/eclipse/rst/aquifer.hpp>
#include <opm/parser/eclipse/Deck/DeckSection.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -46,7 +48,6 @@
#include <opm/parser/eclipse/Units/Dimension.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
namespace Opm {
@@ -266,6 +267,11 @@ AquiferConfig load_aquifers(const Deck& deck, const TableManager& tables, NNC& i
return this->tracer_config;
}
void EclipseState::loadRestartAquifers(const RestartIO::RstAquifer& aquifers) {
if (aquifers.hasAnalyticAquifers())
this->aquifer_config.loadFromRestart(aquifers, this->m_tables);
}
void EclipseState::initTransMult() {
const auto& fp = this->field_props;
if (fp.has_double("MULTX")) this->m_transMult.applyMULT(fp.get_global_double("MULTX") , FaceDir::XPlus);