mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Refactor Aquifer datastructures to follow opm-common
- Replace use of Base:: with this-> in derived classes - Add AquiferInterface::size() utility functions - Remove AquiferInterface::cell_idx_ member
This commit is contained in:
parent
76682497e8
commit
3a7f5799af
@ -49,19 +49,19 @@ public:
|
|||||||
|
|
||||||
using Base::waterCompIdx;
|
using Base::waterCompIdx;
|
||||||
using Base::waterPhaseIdx;
|
using Base::waterPhaseIdx;
|
||||||
AquiferCarterTracy(const Aquancon::AquanconOutput& connection,
|
AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
|
||||||
const std::unordered_map<int, int>& cartesian_to_compressed,
|
const std::unordered_map<int, int>& cartesian_to_compressed,
|
||||||
const Simulator& ebosSimulator,
|
const Simulator& ebosSimulator,
|
||||||
const AquiferCT::AQUCT_data& aquct_data)
|
const AquiferCT::AQUCT_data& aquct_data)
|
||||||
: Base(connection, cartesian_to_compressed, ebosSimulator)
|
: Base(aquct_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
|
||||||
, aquct_data_(aquct_data)
|
, aquct_data_(aquct_data)
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
void endTimeStep() override
|
void endTimeStep() override
|
||||||
{
|
{
|
||||||
for (const auto& Qai : Base::Qai_) {
|
for (const auto& q : this->Qai_) {
|
||||||
Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize();
|
this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -75,17 +75,16 @@ protected:
|
|||||||
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
|
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
|
||||||
inline void initializeConnections() override
|
inline void initializeConnections() override
|
||||||
{
|
{
|
||||||
const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
|
const auto& eclState = this->ebos_simulator_.vanguard().eclState();
|
||||||
const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
|
const auto& ugrid = this->ebos_simulator_.vanguard().grid();
|
||||||
const auto& grid = eclState.getInputGrid();
|
const auto& grid = eclState.getInputGrid();
|
||||||
|
|
||||||
Base::cell_idx_ = this->connection_.global_index;
|
|
||||||
auto globalCellIdx = ugrid.globalCell();
|
auto globalCellIdx = ugrid.globalCell();
|
||||||
|
|
||||||
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
||||||
Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0);
|
this->cell_depth_.resize(this->size(), aquct_data_.d0);
|
||||||
Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
|
this->alphai_.resize(this->size(), 1.0);
|
||||||
Base::faceArea_connected_.resize(Base::cell_idx_.size(), 0.0);
|
this->faceArea_connected_.resize(this->size(), 0.0);
|
||||||
|
|
||||||
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
|
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
|
||||||
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
|
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
|
||||||
@ -95,10 +94,10 @@ protected:
|
|||||||
|
|
||||||
// denom_face_areas is the sum of the areas connected to an aquifer
|
// denom_face_areas is the sum of the areas connected to an aquifer
|
||||||
Scalar denom_face_areas = 0.;
|
Scalar denom_face_areas = 0.;
|
||||||
Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
|
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
|
||||||
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
|
for (size_t idx = 0; idx < this->size(); ++idx) {
|
||||||
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
|
const int cell_index = this->cartesian_to_compressed_.at(this->connections_[idx].global_index);
|
||||||
Base::cellToConnectionIdx_[cell_index] = idx;
|
this->cellToConnectionIdx_[cell_index] = idx;
|
||||||
|
|
||||||
const auto cellFacesRange = cell2Faces[cell_index];
|
const auto cellFacesRange = cell2Faces[cell_index];
|
||||||
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
|
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
|
||||||
@ -132,21 +131,21 @@ protected:
|
|||||||
"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
|
"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
|
||||||
}
|
}
|
||||||
|
|
||||||
if (faceDirection == this->connection_.reservoir_face_dir.at(idx)) {
|
if (faceDirection == this->connections_[idx].face_dir) {
|
||||||
Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx);
|
this->faceArea_connected_.at(idx) = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
|
||||||
denom_face_areas += (this->connection_.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx));
|
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
|
auto cellCenter = grid.getCellCenter(this->connections_[idx].global_index);
|
||||||
Base::cell_depth_.at(idx) = cellCenter[2];
|
this->cell_depth_.at(idx) = cellCenter[2];
|
||||||
}
|
}
|
||||||
|
|
||||||
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
||||||
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
|
for (size_t idx = 0; idx < this->size(); ++idx) {
|
||||||
Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)
|
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
|
||||||
? // Prevent no connection NaNs due to division by zero
|
? // Prevent no connection NaNs due to division by zero
|
||||||
0.
|
0.
|
||||||
: (this->connection_.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas;
|
: (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx)) / denom_face_areas;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -165,22 +164,22 @@ protected:
|
|||||||
|
|
||||||
inline Scalar dpai(int idx)
|
inline Scalar dpai(int idx)
|
||||||
{
|
{
|
||||||
Scalar dp = Base::pa0_
|
Scalar dp = this->pa0_
|
||||||
+ Base::rhow_.at(idx).value() * Base::gravity_() * (Base::cell_depth_.at(idx) - aquct_data_.d0)
|
+ this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - aquct_data_.d0)
|
||||||
- Base::pressure_previous_.at(idx);
|
- this->pressure_previous_.at(idx);
|
||||||
return dp;
|
return dp;
|
||||||
}
|
}
|
||||||
|
|
||||||
// This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
|
// This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
|
||||||
inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator)
|
inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator)
|
||||||
{
|
{
|
||||||
const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_;
|
const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_;
|
||||||
const Scalar td = simulator.time() / Base::Tc_;
|
const Scalar td = simulator.time() / this->Tc_;
|
||||||
Scalar PItdprime = 0.;
|
Scalar PItdprime = 0.;
|
||||||
Scalar PItd = 0.;
|
Scalar PItd = 0.;
|
||||||
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
|
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
|
||||||
a = 1.0 / Base::Tc_ * ((beta_ * dpai(idx)) - (Base::W_flux_.value() * PItdprime)) / (PItd - td * PItdprime);
|
a = 1.0 / this->Tc_ * ((beta_ * dpai(idx)) - (this->W_flux_.value() * PItdprime)) / (PItd - td * PItdprime);
|
||||||
b = beta_ / (Base::Tc_ * (PItd - td * PItdprime));
|
b = beta_ / (this->Tc_ * (PItd - td * PItdprime));
|
||||||
}
|
}
|
||||||
|
|
||||||
// This function implements Eq 5.7 of the EclipseTechnicalDescription
|
// This function implements Eq 5.7 of the EclipseTechnicalDescription
|
||||||
@ -188,8 +187,8 @@ protected:
|
|||||||
{
|
{
|
||||||
Scalar a, b;
|
Scalar a, b;
|
||||||
calculateEqnConstants(a, b, idx, simulator);
|
calculateEqnConstants(a, b, idx, simulator);
|
||||||
Base::Qai_.at(idx)
|
this->Qai_.at(idx)
|
||||||
= Base::alphai_.at(idx) * (a - b * (Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx)));
|
= this->alphai_.at(idx) * (a - b * (this->pressure_current_.at(idx) - this->pressure_previous_.at(idx)));
|
||||||
}
|
}
|
||||||
|
|
||||||
inline void calculateAquiferConstants() override
|
inline void calculateAquiferConstants() override
|
||||||
@ -198,7 +197,7 @@ protected:
|
|||||||
beta_ = aquct_data_.c2 * aquct_data_.h * aquct_data_.theta * aquct_data_.phi_aq * aquct_data_.C_t
|
beta_ = aquct_data_.c2 * aquct_data_.h * aquct_data_.theta * aquct_data_.phi_aq * aquct_data_.C_t
|
||||||
* aquct_data_.r_o * aquct_data_.r_o;
|
* aquct_data_.r_o * aquct_data_.r_o;
|
||||||
// We calculate the time constant
|
// We calculate the time constant
|
||||||
Base::Tc_ = mu_w_ * aquct_data_.phi_aq * aquct_data_.C_t * aquct_data_.r_o * aquct_data_.r_o
|
this->Tc_ = mu_w_ * aquct_data_.phi_aq * aquct_data_.C_t * aquct_data_.r_o * aquct_data_.r_o
|
||||||
/ (aquct_data_.k_a * aquct_data_.c1);
|
/ (aquct_data_.k_a * aquct_data_.c1);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -206,17 +205,17 @@ protected:
|
|||||||
{
|
{
|
||||||
|
|
||||||
int pvttableIdx = aquct_data_.pvttableID - 1;
|
int pvttableIdx = aquct_data_.pvttableID - 1;
|
||||||
Base::rhow_.resize(Base::cell_idx_.size(), 0.);
|
this->rhow_.resize(this->size(), 0.);
|
||||||
if (!aquct_data_.p0.first) {
|
if (!aquct_data_.p0.first) {
|
||||||
Base::pa0_ = calculateReservoirEquilibrium();
|
this->pa0_ = calculateReservoirEquilibrium();
|
||||||
} else {
|
} else {
|
||||||
Base::pa0_ = aquct_data_.p0.second;
|
this->pa0_ = aquct_data_.p0.second;
|
||||||
}
|
}
|
||||||
|
|
||||||
// use the thermodynamic state of the first active cell as a
|
// use the thermodynamic state of the first active cell as a
|
||||||
// reference. there might be better ways to do this...
|
// reference. there might be better ways to do this...
|
||||||
ElementContext elemCtx(Base::ebos_simulator_);
|
ElementContext elemCtx(this->ebos_simulator_);
|
||||||
auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>();
|
auto elemIt = this->ebos_simulator_.gridView().template begin</*codim=*/0>();
|
||||||
elemCtx.updatePrimaryStencil(*elemIt);
|
elemCtx.updatePrimaryStencil(*elemIt);
|
||||||
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
|
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
|
||||||
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
|
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
|
||||||
@ -227,7 +226,7 @@ protected:
|
|||||||
fs_aquifer.assign(iq0.fluidState());
|
fs_aquifer.assign(iq0.fluidState());
|
||||||
Eval temperature_aq, pa0_mean;
|
Eval temperature_aq, pa0_mean;
|
||||||
temperature_aq = fs_aquifer.temperature(0);
|
temperature_aq = fs_aquifer.temperature(0);
|
||||||
pa0_mean = Base::pa0_;
|
pa0_mean = this->pa0_;
|
||||||
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
|
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
|
||||||
mu_w_ = mu_w_aquifer.value();
|
mu_w_ = mu_w_aquifer.value();
|
||||||
}
|
}
|
||||||
@ -240,8 +239,8 @@ protected:
|
|||||||
std::vector<Scalar> pw_aquifer;
|
std::vector<Scalar> pw_aquifer;
|
||||||
Scalar water_pressure_reservoir;
|
Scalar water_pressure_reservoir;
|
||||||
|
|
||||||
ElementContext elemCtx(Base::ebos_simulator_);
|
ElementContext elemCtx(this->ebos_simulator_);
|
||||||
const auto& gridView = Base::ebos_simulator_.gridView();
|
const auto& gridView = this->ebos_simulator_.gridView();
|
||||||
auto elemIt = gridView.template begin</*codim=*/0>();
|
auto elemIt = gridView.template begin</*codim=*/0>();
|
||||||
const auto& elemEndIt = gridView.template end</*codim=*/0>();
|
const auto& elemEndIt = gridView.template end</*codim=*/0>();
|
||||||
for (; elemIt != elemEndIt; ++elemIt) {
|
for (; elemIt != elemEndIt; ++elemIt) {
|
||||||
@ -249,7 +248,7 @@ protected:
|
|||||||
elemCtx.updatePrimaryStencil(elem);
|
elemCtx.updatePrimaryStencil(elem);
|
||||||
|
|
||||||
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
|
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
|
||||||
int idx = Base::cellToConnectionIdx_[cellIdx];
|
int idx = this->cellToConnectionIdx_[cellIdx];
|
||||||
if (idx < 0)
|
if (idx < 0)
|
||||||
continue;
|
continue;
|
||||||
|
|
||||||
@ -258,11 +257,11 @@ protected:
|
|||||||
const auto& fs = iq0.fluidState();
|
const auto& fs = iq0.fluidState();
|
||||||
|
|
||||||
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
|
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
|
||||||
Base::rhow_[idx] = fs.density(waterPhaseIdx);
|
this->rhow_[idx] = fs.density(waterPhaseIdx);
|
||||||
pw_aquifer.push_back(
|
pw_aquifer.push_back(
|
||||||
(water_pressure_reservoir
|
(water_pressure_reservoir
|
||||||
- Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aquct_data_.d0))
|
- this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - aquct_data_.d0))
|
||||||
* Base::alphai_[idx]);
|
* this->alphai_[idx]);
|
||||||
}
|
}
|
||||||
|
|
||||||
// We take the average of the calculated equilibrium pressures.
|
// We take the average of the calculated equilibrium pressures.
|
||||||
|
@ -51,19 +51,19 @@ public:
|
|||||||
using Base::waterCompIdx;
|
using Base::waterCompIdx;
|
||||||
using Base::waterPhaseIdx;
|
using Base::waterPhaseIdx;
|
||||||
|
|
||||||
AquiferFetkovich(const Aquancon::AquanconOutput& connection,
|
AquiferFetkovich(const std::vector<Aquancon::AquancCell>& connections,
|
||||||
const std::unordered_map<int, int>& cartesian_to_compressed,
|
const std::unordered_map<int, int>& cartesian_to_compressed,
|
||||||
const Simulator& ebosSimulator,
|
const Simulator& ebosSimulator,
|
||||||
const Aquifetp::AQUFETP_data& aqufetp_data)
|
const Aquifetp::AQUFETP_data& aqufetp_data)
|
||||||
: Base(connection, cartesian_to_compressed, ebosSimulator)
|
: Base(aqufetp_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
|
||||||
, aqufetp_data_(aqufetp_data)
|
, aqufetp_data_(aqufetp_data)
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
void endTimeStep() override
|
void endTimeStep() override
|
||||||
{
|
{
|
||||||
for (const auto& Qai : Base::Qai_) {
|
for (const auto& q : this->Qai_) {
|
||||||
Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize();
|
this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
|
||||||
aquifer_pressure_ = aquiferPressure();
|
aquifer_pressure_ = aquiferPressure();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -76,18 +76,17 @@ protected:
|
|||||||
|
|
||||||
inline void initializeConnections() override
|
inline void initializeConnections() override
|
||||||
{
|
{
|
||||||
const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
|
const auto& eclState = this->ebos_simulator_.vanguard().eclState();
|
||||||
const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
|
const auto& ugrid = this->ebos_simulator_.vanguard().grid();
|
||||||
const auto& grid = eclState.getInputGrid();
|
const auto& grid = eclState.getInputGrid();
|
||||||
|
|
||||||
Base::cell_idx_ = this->connection_.global_index;
|
|
||||||
auto globalCellIdx = ugrid.globalCell();
|
auto globalCellIdx = ugrid.globalCell();
|
||||||
|
|
||||||
|
|
||||||
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
||||||
Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0);
|
this->cell_depth_.resize(this->size(), aqufetp_data_.d0);
|
||||||
Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
|
this->alphai_.resize(this->size(), 1.0);
|
||||||
Base::faceArea_connected_.resize(Base::cell_idx_.size(), 0.0);
|
this->faceArea_connected_.resize(this->size(), 0.0);
|
||||||
|
|
||||||
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
|
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
|
||||||
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
|
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
|
||||||
@ -97,14 +96,16 @@ protected:
|
|||||||
|
|
||||||
// denom_face_areas is the sum of the areas connected to an aquifer
|
// denom_face_areas is the sum of the areas connected to an aquifer
|
||||||
Scalar denom_face_areas = 0.;
|
Scalar denom_face_areas = 0.;
|
||||||
Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
|
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
|
||||||
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
|
for (size_t idx = 0; idx < this->size(); ++idx) {
|
||||||
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
|
const auto global_index = this->connections_[idx].global_index;
|
||||||
Base::cellToConnectionIdx_[cell_index] = idx;
|
const int cell_index = this->cartesian_to_compressed_.at(global_index);
|
||||||
const auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
|
|
||||||
Base::cell_depth_.at(idx) = cellCenter[2];
|
|
||||||
|
|
||||||
if (!this->connection_.influx_coeff[idx]) { // influx_coeff is defaulted
|
this->cellToConnectionIdx_[cell_index] = idx;
|
||||||
|
const auto cellCenter = grid.getCellCenter(global_index);
|
||||||
|
this->cell_depth_.at(idx) = cellCenter[2];
|
||||||
|
|
||||||
|
if (!this->connections_[idx].influx_coeff.first) { // influx_coeff is defaulted
|
||||||
const auto cellFacesRange = cell2Faces[cell_index];
|
const auto cellFacesRange = cell2Faces[cell_index];
|
||||||
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
|
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
|
||||||
// The index of the face in the compressed grid
|
// The index of the face in the compressed grid
|
||||||
@ -137,24 +138,23 @@ protected:
|
|||||||
"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
|
"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
|
||||||
}
|
}
|
||||||
|
|
||||||
if (faceDirection == this->connection_.reservoir_face_dir.at(idx)) {
|
if (faceDirection == this->connections_[idx].face_dir) {
|
||||||
Base::faceArea_connected_.at(idx)
|
this->faceArea_connected_[idx] = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
|
||||||
= Base::getFaceArea(faceCells, ugrid, faceIdx, idx);
|
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
Base::faceArea_connected_.at(idx) = *this->connection_.influx_coeff[idx];
|
this->faceArea_connected_.at(idx) = this->connections_[idx].influx_coeff.second;
|
||||||
}
|
}
|
||||||
denom_face_areas += (this->connection_.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx));
|
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
|
||||||
}
|
}
|
||||||
|
|
||||||
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
||||||
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
|
for (size_t idx = 0; idx < this->size(); ++idx) {
|
||||||
Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)
|
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
|
||||||
? // Prevent no connection NaNs due to division by zero
|
? // Prevent no connection NaNs due to division by zero
|
||||||
0.
|
0.
|
||||||
: (this->connection_.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas;
|
: (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx)) / denom_face_areas;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -170,45 +170,45 @@ protected:
|
|||||||
|
|
||||||
inline Eval dpai(int idx)
|
inline Eval dpai(int idx)
|
||||||
{
|
{
|
||||||
const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx)
|
const Eval dp = aquifer_pressure_ - this->pressure_current_.at(idx)
|
||||||
+ Base::rhow_[idx] * Base::gravity_() * (Base::cell_depth_[idx] - aqufetp_data_.d0);
|
+ this->rhow_[idx] * this->gravity_() * (this->cell_depth_[idx] - aqufetp_data_.d0);
|
||||||
return dp;
|
return dp;
|
||||||
}
|
}
|
||||||
|
|
||||||
// This function implements Eq 5.12 of the EclipseTechnicalDescription
|
// This function implements Eq 5.12 of the EclipseTechnicalDescription
|
||||||
inline Scalar aquiferPressure()
|
inline Scalar aquiferPressure()
|
||||||
{
|
{
|
||||||
Scalar Flux = Base::W_flux_.value();
|
Scalar Flux = this->W_flux_.value();
|
||||||
Scalar pa_ = Base::pa0_ - Flux / (aqufetp_data_.C_t * aqufetp_data_.V0);
|
Scalar pa_ = this->pa0_ - Flux / (aqufetp_data_.C_t * aqufetp_data_.V0);
|
||||||
return pa_;
|
return pa_;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline void calculateAquiferConstants() override
|
inline void calculateAquiferConstants() override
|
||||||
{
|
{
|
||||||
Base::Tc_ = (aqufetp_data_.C_t * aqufetp_data_.V0) / aqufetp_data_.J;
|
this->Tc_ = (aqufetp_data_.C_t * aqufetp_data_.V0) / aqufetp_data_.J;
|
||||||
}
|
}
|
||||||
// This function implements Eq 5.14 of the EclipseTechnicalDescription
|
// This function implements Eq 5.14 of the EclipseTechnicalDescription
|
||||||
inline void calculateInflowRate(int idx, const Simulator& simulator) override
|
inline void calculateInflowRate(int idx, const Simulator& simulator) override
|
||||||
{
|
{
|
||||||
const Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_;
|
const Scalar td_Tc_ = simulator.timeStepSize() / this->Tc_;
|
||||||
const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_;
|
const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_;
|
||||||
Base::Qai_.at(idx) = Base::alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef;
|
this->Qai_.at(idx) = this->alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline void calculateAquiferCondition() override
|
inline void calculateAquiferCondition() override
|
||||||
{
|
{
|
||||||
Base::rhow_.resize(Base::cell_idx_.size(), 0.);
|
this->rhow_.resize(this->size(), 0.);
|
||||||
|
|
||||||
if (this->solution_set_from_restart_) {
|
if (this->solution_set_from_restart_) {
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
if (!aqufetp_data_.p0.first) {
|
if (!aqufetp_data_.p0.first) {
|
||||||
Base::pa0_ = calculateReservoirEquilibrium();
|
this->pa0_ = calculateReservoirEquilibrium();
|
||||||
} else {
|
} else {
|
||||||
Base::pa0_ = aqufetp_data_.p0.second;
|
this->pa0_ = aqufetp_data_.p0.second;
|
||||||
}
|
}
|
||||||
aquifer_pressure_ = Base::pa0_;
|
aquifer_pressure_ = this->pa0_;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline Scalar calculateReservoirEquilibrium() override
|
inline Scalar calculateReservoirEquilibrium() override
|
||||||
@ -217,15 +217,15 @@ protected:
|
|||||||
std::vector<Scalar> pw_aquifer;
|
std::vector<Scalar> pw_aquifer;
|
||||||
Scalar water_pressure_reservoir;
|
Scalar water_pressure_reservoir;
|
||||||
|
|
||||||
ElementContext elemCtx(Base::ebos_simulator_);
|
ElementContext elemCtx(this->ebos_simulator_);
|
||||||
const auto& gridView = Base::ebos_simulator_.gridView();
|
const auto& gridView = this->ebos_simulator_.gridView();
|
||||||
auto elemIt = gridView.template begin</*codim=*/0>();
|
auto elemIt = gridView.template begin</*codim=*/0>();
|
||||||
const auto& elemEndIt = gridView.template end</*codim=*/0>();
|
const auto& elemEndIt = gridView.template end</*codim=*/0>();
|
||||||
for (; elemIt != elemEndIt; ++elemIt) {
|
for (; elemIt != elemEndIt; ++elemIt) {
|
||||||
const auto& elem = *elemIt;
|
const auto& elem = *elemIt;
|
||||||
elemCtx.updatePrimaryStencil(elem);
|
elemCtx.updatePrimaryStencil(elem);
|
||||||
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
|
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
|
||||||
int idx = Base::cellToConnectionIdx_[cellIdx];
|
int idx = this->cellToConnectionIdx_[cellIdx];
|
||||||
if (idx < 0)
|
if (idx < 0)
|
||||||
continue;
|
continue;
|
||||||
|
|
||||||
@ -234,11 +234,11 @@ protected:
|
|||||||
const auto& fs = iq0.fluidState();
|
const auto& fs = iq0.fluidState();
|
||||||
|
|
||||||
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
|
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
|
||||||
Base::rhow_[idx] = fs.density(waterPhaseIdx);
|
this->rhow_[idx] = fs.density(waterPhaseIdx);
|
||||||
pw_aquifer.push_back(
|
pw_aquifer.push_back(
|
||||||
(water_pressure_reservoir
|
(water_pressure_reservoir
|
||||||
- Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aqufetp_data_.d0))
|
- this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - aqufetp_data_.d0))
|
||||||
* Base::alphai_[idx]);
|
* this->alphai_[idx]);
|
||||||
}
|
}
|
||||||
|
|
||||||
// We take the average of the calculated equilibrium pressures.
|
// We take the average of the calculated equilibrium pressures.
|
||||||
|
@ -73,16 +73,15 @@ public:
|
|||||||
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
|
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
|
||||||
|
|
||||||
// Constructor
|
// Constructor
|
||||||
AquiferInterface(const Aquancon::AquanconOutput& connection,
|
AquiferInterface(int aqID,
|
||||||
|
const std::vector<Aquancon::AquancCell>& connections,
|
||||||
const std::unordered_map<int, int>& cartesian_to_compressed,
|
const std::unordered_map<int, int>& cartesian_to_compressed,
|
||||||
const Simulator& ebosSimulator)
|
const Simulator& ebosSimulator)
|
||||||
: connection_(connection)
|
: aquiferID(aqID)
|
||||||
|
, connections_(connections)
|
||||||
, ebos_simulator_(ebosSimulator)
|
, ebos_simulator_(ebosSimulator)
|
||||||
, cartesian_to_compressed_(cartesian_to_compressed)
|
, cartesian_to_compressed_(cartesian_to_compressed)
|
||||||
{
|
{
|
||||||
assert(this->connection_.influx_coeff.size() == this->connection_.global_index.size());
|
|
||||||
assert(this->connection_.influx_coeff.size() == this->connection_.influx_multiplier.size());
|
|
||||||
assert(this->connection_.influx_multiplier.size() == this->connection_.reservoir_face_dir.size());
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Deconstructor
|
// Deconstructor
|
||||||
@ -94,16 +93,13 @@ public:
|
|||||||
{
|
{
|
||||||
auto xaqPos
|
auto xaqPos
|
||||||
= std::find_if(aquiferSoln.begin(), aquiferSoln.end(), [this](const data::AquiferData& xaq) -> bool {
|
= std::find_if(aquiferSoln.begin(), aquiferSoln.end(), [this](const data::AquiferData& xaq) -> bool {
|
||||||
return xaq.aquiferID == this->connection_.aquiferID;
|
return xaq.aquiferID == this->aquiferID;
|
||||||
});
|
});
|
||||||
|
|
||||||
if (xaqPos == aquiferSoln.end()) {
|
if (xaqPos == aquiferSoln.end())
|
||||||
// No restart value applies to this aquifer. Nothing to do.
|
|
||||||
return;
|
return;
|
||||||
}
|
|
||||||
|
|
||||||
this->assignRestartData(*xaqPos);
|
this->assignRestartData(*xaqPos);
|
||||||
|
|
||||||
this->W_flux_ = xaqPos->volume;
|
this->W_flux_ = xaqPos->volume;
|
||||||
this->pa0_ = xaqPos->initPressure;
|
this->pa0_ = xaqPos->initPressure;
|
||||||
this->solution_set_from_restart_ = true;
|
this->solution_set_from_restart_ = true;
|
||||||
@ -155,6 +151,12 @@ public:
|
|||||||
+= Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
|
+= Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
std::size_t size() const {
|
||||||
|
return this->connections_.size();
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
inline Scalar gravity_() const
|
inline Scalar gravity_() const
|
||||||
{
|
{
|
||||||
@ -174,9 +176,9 @@ protected:
|
|||||||
calculateAquiferCondition();
|
calculateAquiferCondition();
|
||||||
calculateAquiferConstants();
|
calculateAquiferConstants();
|
||||||
|
|
||||||
pressure_previous_.resize(cell_idx_.size(), 0.);
|
pressure_previous_.resize(this->connections_.size(), 0.);
|
||||||
pressure_current_.resize(cell_idx_.size(), 0.);
|
pressure_current_.resize(this->connections_.size(), 0.);
|
||||||
Qai_.resize(cell_idx_.size(), 0.0);
|
Qai_.resize(this->connections_.size(), 0.0);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline void
|
inline void
|
||||||
@ -214,7 +216,7 @@ protected:
|
|||||||
const auto cellNeighbour1 = faceCells(faceIdx, 1);
|
const auto cellNeighbour1 = faceCells(faceIdx, 1);
|
||||||
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
|
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
|
||||||
const auto calculatedFaceArea
|
const auto calculatedFaceArea
|
||||||
= (!this->connection_.influx_coeff.at(idx)) ? defaultFaceArea : *(this->connection_.influx_coeff.at(idx));
|
= (!this->connections_[idx].influx_coeff.first) ? defaultFaceArea : this->connections_[idx].influx_coeff.second;
|
||||||
faceArea = (cellNeighbour0 * cellNeighbour1 > 0) ? 0. : calculatedFaceArea;
|
faceArea = (cellNeighbour0 * cellNeighbour1 > 0) ? 0. : calculatedFaceArea;
|
||||||
if (cellNeighbour1 == 0) {
|
if (cellNeighbour1 == 0) {
|
||||||
faceArea = (cellNeighbour0 < 0) ? faceArea : 0.;
|
faceArea = (cellNeighbour0 < 0) ? faceArea : 0.;
|
||||||
@ -226,12 +228,12 @@ protected:
|
|||||||
|
|
||||||
virtual void endTimeStep() = 0;
|
virtual void endTimeStep() = 0;
|
||||||
|
|
||||||
const Aquancon::AquanconOutput connection_;
|
const int aquiferID;
|
||||||
|
const std::vector<Aquancon::AquancCell> connections_;
|
||||||
const Simulator& ebos_simulator_;
|
const Simulator& ebos_simulator_;
|
||||||
const std::unordered_map<int, int> cartesian_to_compressed_;
|
const std::unordered_map<int, int> cartesian_to_compressed_;
|
||||||
|
|
||||||
// Grid variables
|
// Grid variables
|
||||||
std::vector<size_t> cell_idx_;
|
|
||||||
std::vector<Scalar> faceArea_connected_;
|
std::vector<Scalar> faceArea_connected_;
|
||||||
std::vector<int> cellToConnectionIdx_;
|
std::vector<int> cellToConnectionIdx_;
|
||||||
// Quantities at each grid id
|
// Quantities at each grid id
|
||||||
|
@ -158,20 +158,15 @@ BlackoilAquiferModel<TypeTag>::init()
|
|||||||
const AquiferCT aquiferct = AquiferCT(eclState.getTableManager(), deck);
|
const AquiferCT aquiferct = AquiferCT(eclState.getTableManager(), deck);
|
||||||
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
||||||
|
|
||||||
std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.data();
|
|
||||||
std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
|
|
||||||
|
|
||||||
assert(aquifersData.size() == aquifer_connection.size());
|
|
||||||
const auto& ugrid = simulator_.vanguard().grid();
|
const auto& ugrid = simulator_.vanguard().grid();
|
||||||
const auto& gridView = simulator_.gridView();
|
const auto& gridView = simulator_.gridView();
|
||||||
const int number_of_cells = gridView.size(0);
|
const int number_of_cells = gridView.size(0);
|
||||||
|
|
||||||
cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
|
cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
|
||||||
|
|
||||||
for (size_t i = 0; i < aquifersData.size(); ++i) {
|
for (const auto& aquifer : aquiferct)
|
||||||
aquifers_CarterTracy.push_back(AquiferCarterTracy<TypeTag>(
|
aquifers_CarterTracy.push_back(AquiferCarterTracy<TypeTag>(aquifer_connect[aquifer.aquiferID], cartesian_to_compressed_, this->simulator_, aquifer));
|
||||||
aquifer_connection[i], cartesian_to_compressed_, this->simulator_, aquifersData[i]));
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
if (comm.rank() == 0)
|
if (comm.rank() == 0)
|
||||||
has = deck.hasKeyword("AQUFETP");
|
has = deck.hasKeyword("AQUFETP");
|
||||||
@ -188,20 +183,15 @@ BlackoilAquiferModel<TypeTag>::init()
|
|||||||
const Aquifetp aquifetp = Aquifetp(deck);
|
const Aquifetp aquifetp = Aquifetp(deck);
|
||||||
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
||||||
|
|
||||||
std::vector<Aquifetp::AQUFETP_data> aquifersData = aquifetp.data();
|
|
||||||
std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
|
|
||||||
|
|
||||||
assert(aquifersData.size() == aquifer_connection.size());
|
|
||||||
const auto& ugrid = simulator_.vanguard().grid();
|
const auto& ugrid = simulator_.vanguard().grid();
|
||||||
const auto& gridView = simulator_.gridView();
|
const auto& gridView = simulator_.gridView();
|
||||||
const int number_of_cells = gridView.size(0);
|
const int number_of_cells = gridView.size(0);
|
||||||
|
|
||||||
cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
|
cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
|
||||||
|
|
||||||
for (size_t i = 0; i < aquifersData.size(); ++i) {
|
for (const auto& aquifer : aquifetp)
|
||||||
aquifers_Fetkovich.push_back(AquiferFetkovich<TypeTag>(
|
aquifers_Fetkovich.push_back(AquiferFetkovich<TypeTag>(aquifer_connect[aquifer.aquiferID], cartesian_to_compressed_, this->simulator_, aquifer));
|
||||||
aquifer_connection[i], cartesian_to_compressed_, this->simulator_, aquifersData[i]));
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
template <typename TypeTag>
|
template <typename TypeTag>
|
||||||
|
Loading…
Reference in New Issue
Block a user