OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2254 from GitPaean/clang-format-aquifer

Browse Source 
re-formatting aquifer files with clang-format
This commit is contained in:
Bård Skaflestad 2019-12-20 16:50:11 +01:00 committed by GitHub
commit b64116a70c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 679 additions and 685 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

421
opm/simulators/aquifers/AquiferCarterTracy.hpp
View File

@ -31,249 +31,250 @@
namespace Opm namespace Opm
{ {
template<typename TypeTag> template <typename TypeTag>
class AquiferCarterTracy: public AquiferInterface<TypeTag> class AquiferCarterTracy : public AquiferInterface<TypeTag>
{
public:
typedef AquiferInterface<TypeTag> Base;
using typename Base::BlackoilIndices;
using typename Base::ElementContext;
using typename Base::Eval;
using typename Base::FluidState;
using typename Base::FluidSystem;
using typename Base::IntensiveQuantities;
using typename Base::RateVector;
using typename Base::Scalar;
using typename Base::Simulator;
using Base::waterCompIdx;
using Base::waterPhaseIdx;
AquiferCarterTracy(const Aquancon::AquanconOutput& connection,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator,
const AquiferCT::AQUCT_data& aquct_data)
: Base(connection, cartesian_to_compressed, ebosSimulator)
, aquct_data_(aquct_data)
{ {
public: }
typedef AquiferInterface<TypeTag> Base;
using typename Base::Simulator; void endTimeStep() override
using typename Base::ElementContext; {
using typename Base::FluidSystem; for (const auto& Qai : Base::Qai_) {
using typename Base::BlackoilIndices; Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize();
using typename Base::RateVector; }
using typename Base::IntensiveQuantities; }
using typename Base::Eval;
using typename Base::Scalar;
using typename Base::FluidState;
using Base::waterCompIdx; protected:
using Base::waterPhaseIdx; // Variables constants
AquiferCarterTracy( const Aquancon::AquanconOutput& connection, const AquiferCT::AQUCT_data aquct_data_;
const std::unordered_map<int, int>& cartesian_to_compressed, Scalar beta_; // Influx constant
const Simulator& ebosSimulator, // TODO: it is possible it should be a AD variable
const AquiferCT::AQUCT_data& aquct_data) Scalar mu_w_; // water viscosity
: Base(connection, cartesian_to_compressed, ebosSimulator)
, aquct_data_(aquct_data)
{}
void endTimeStep() override // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
{ inline void initializeConnections(const Aquancon::AquanconOutput& connection) override
for (const auto& Qai: Base::Qai_) { {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize(); const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
} const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
} const auto& grid = eclState.getInputGrid();
protected: Base::cell_idx_ = connection.global_index;
// Variables constants auto globalCellIdx = ugrid.globalCell();
const AquiferCT::AQUCT_data aquct_data_;
Scalar beta_; // Influx constant
// TODO: it is possible it should be a AD variable
Scalar mu_w_; // water viscosity
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer assert(Base::cell_idx_ == connection.global_index);
inline void initializeConnections(const Aquancon::AquanconOutput& connection) override assert((Base::cell_idx_.size() <= connection.influx_coeff.size()));
{ assert((connection.influx_coeff.size() == connection.influx_multiplier.size()));
const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); assert((connection.influx_multiplier.size() == connection.reservoir_face_dir.size()));
const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
Base::cell_idx_ = connection.global_index; // We hack the cell depth values for now. We can actually get it from elementcontext pos
auto globalCellIdx = ugrid.globalCell(); Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0);
Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
Base::faceArea_connected_.resize(Base::cell_idx_.size(), 0.0);
assert( Base::cell_idx_ == connection.global_index); auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
assert( (Base::cell_idx_.size() <= connection.influx_coeff.size()) ); auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) );
assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) );
// We hack the cell depth values for now. We can actually get it from elementcontext pos // Translate the C face tag into the enum used by opm-parser's TransMult class
Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0); Opm::FaceDir::DirEnum faceDirection;
Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); // denom_face_areas is the sum of the areas connected to an aquifer
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); Scalar denom_face_areas = 0.;
Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
Base::cellToConnectionIdx_[cell_index] = idx;
// Translate the C face tag into the enum used by opm-parser's TransMult class const auto cellFacesRange = cell2Faces[cell_index];
Opm::FaceDir::DirEnum faceDirection; for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
// The index of the face in the compressed grid
const int faceIdx = *cellFaceIter;
// denom_face_areas is the sum of the areas connected to an aquifer // the logically-Cartesian direction of the face
Scalar denom_face_areas = 0.; const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx)
{
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
Base::cellToConnectionIdx_[cell_index] = idx;
const auto cellFacesRange = cell2Faces[cell_index]; switch (faceTag) {
for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) case 0:
{ faceDirection = Opm::FaceDir::XMinus;
// The index of the face in the compressed grid break;
const int faceIdx = *cellFaceIter; case 1:
faceDirection = Opm::FaceDir::XPlus;
// the logically-Cartesian direction of the face break;
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); case 2:
faceDirection = Opm::FaceDir::YMinus;
switch(faceTag) break;
{ case 3:
case 0: faceDirection = Opm::FaceDir::XMinus; faceDirection = Opm::FaceDir::YPlus;
break; break;
case 1: faceDirection = Opm::FaceDir::XPlus; case 4:
break; faceDirection = Opm::FaceDir::ZMinus;
case 2: faceDirection = Opm::FaceDir::YMinus; break;
break; case 5:
case 3: faceDirection = Opm::FaceDir::YPlus; faceDirection = Opm::FaceDir::ZPlus;
break; break;
case 4: faceDirection = Opm::FaceDir::ZMinus; default:
break; OPM_THROW(Opm::NumericalIssue,
case 5: faceDirection = Opm::FaceDir::ZPlus; "Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
break;
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
}
if (faceDirection == connection.reservoir_face_dir.at(idx))
{
Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection);
denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) );
}
}
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::cell_depth_.at(idx) = cellCenter[2];
} }
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon()); if (faceDirection == connection.reservoir_face_dir.at(idx)) {
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection);
{ denom_face_areas += (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx));
Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
0.
: ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas;
} }
} }
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::cell_depth_.at(idx) = cellCenter[2];
}
void assignRestartData(const data::AquiferData& /* xaq */) override const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
{ for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
throw std::runtime_error { Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)
"Restart-based initialization not currently supported " ? // Prevent no connection NaNs due to division by zero
"for Carter-Tracey analytic aquifers" 0.
}; : (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas;
} }
}
inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) void assignRestartData(const data::AquiferData& /* xaq */) override
{ {
// We use the opm-common numeric linear interpolator throw std::runtime_error {"Restart-based initialization not currently supported "
pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); "for Carter-Tracey analytic aquifers"};
pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); }
}
inline Scalar dpai(int idx) inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
{ {
Scalar dp = Base::pa0_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aquct_data_.d0) - Base::pressure_previous_.at(idx); // We use the opm-common numeric linear interpolator
return dp; pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td);
} pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td);
}
// This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription inline Scalar dpai(int idx)
inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) {
{ Scalar dp = Base::pa0_
const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_; + Base::rhow_.at(idx).value() * Base::gravity_() * (Base::cell_depth_.at(idx) - aquct_data_.d0)
const Scalar td = simulator.time() / Base::Tc_; - Base::pressure_previous_.at(idx);
Scalar PItdprime = 0.; return dp;
Scalar PItd = 0.; }
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
a = 1.0/Base::Tc_ * ( (beta_ * dpai(idx)) - (Base::W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime );
b = beta_ / (Base::Tc_ * ( PItd - td*PItdprime));
}
// This function implements Eq 5.7 of the EclipseTechnicalDescription // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator) override inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator)
{ {
Scalar a, b; const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_;
calculateEqnConstants(a,b,idx,simulator); const Scalar td = simulator.time() / Base::Tc_;
Base::Qai_.at(idx) = Base::alphai_.at(idx)*( a - b * ( Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx) ) ); Scalar PItdprime = 0.;
} Scalar PItd = 0.;
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
a = 1.0 / Base::Tc_ * ((beta_ * dpai(idx)) - (Base::W_flux_.value() * PItdprime)) / (PItd - td * PItdprime);
b = beta_ / (Base::Tc_ * (PItd - td * PItdprime));
}
inline void calculateAquiferConstants() override // This function implements Eq 5.7 of the EclipseTechnicalDescription
{ inline void calculateInflowRate(int idx, const Simulator& simulator) override
// We calculate the influx constant {
beta_ = aquct_data_.c2 * aquct_data_.h Scalar a, b;
* aquct_data_.theta * aquct_data_.phi_aq calculateEqnConstants(a, b, idx, simulator);
* aquct_data_.C_t Base::Qai_.at(idx)
* aquct_data_.r_o * aquct_data_.r_o; = Base::alphai_.at(idx) * (a - b * (Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx)));
// 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
/ ( aquct_data_.k_a * aquct_data_.c1 );
}
inline void calculateAquiferCondition() override inline void calculateAquiferConstants() override
{ {
// We calculate the influx constant
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;
// 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
/ (aquct_data_.k_a * aquct_data_.c1);
}
int pvttableIdx = aquct_data_.pvttableID - 1; inline void calculateAquiferCondition() override
Base::rhow_.resize(Base::cell_idx_.size(),0.); {
if (!aquct_data_.p0)
{
Base::pa0_ = calculateReservoirEquilibrium();
}
else
{
Base::pa0_ = *(aquct_data_.p0);
}
// use the thermodynamic state of the first active cell as a int pvttableIdx = aquct_data_.pvttableID - 1;
// reference. there might be better ways to do this... Base::rhow_.resize(Base::cell_idx_.size(), 0.);
ElementContext elemCtx(Base::ebos_simulator_); if (!aquct_data_.p0) {
auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>(); Base::pa0_ = calculateReservoirEquilibrium();
elemCtx.updatePrimaryStencil(*elemIt); } else {
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); Base::pa0_ = *(aquct_data_.p0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); }
// Initialize a FluidState object first
FluidState fs_aquifer;
// We use the temperature of the first cell connected to the aquifer
// Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
fs_aquifer.assign( iq0.fluidState() );
Eval temperature_aq, pa0_mean;
temperature_aq = fs_aquifer.temperature(0);
pa0_mean = Base::pa0_;
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
mu_w_ = mu_w_aquifer.value();
} // use the thermodynamic state of the first active cell as a
// reference. there might be better ways to do this...
ElementContext elemCtx(Base::ebos_simulator_);
auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>();
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
// Initialize a FluidState object first
FluidState fs_aquifer;
// We use the temperature of the first cell connected to the aquifer
// Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
fs_aquifer.assign(iq0.fluidState());
Eval temperature_aq, pa0_mean;
temperature_aq = fs_aquifer.temperature(0);
pa0_mean = Base::pa0_;
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
mu_w_ = mu_w_aquifer.value();
}
// This function is for calculating the aquifer properties from equilibrium state with the reservoir // This function is for calculating the aquifer properties from equilibrium state with the reservoir
// TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models // TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models
inline Scalar calculateReservoirEquilibrium() override inline Scalar calculateReservoirEquilibrium() override
{ {
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
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(Base::ebos_simulator_);
const auto& gridView = Base::ebos_simulator_.gridView(); const auto& gridView = Base::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 = Base::cellToConnectionIdx_[cellIdx];
if (idx < 0) if (idx < 0)
continue; continue;
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);
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); Base::rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back( (water_pressure_reservoir - Base::rhow_[idx].value()*Base::gravity_()*(Base::cell_depth_[idx] - aquct_data_.d0))*Base::alphai_[idx] ); pw_aquifer.push_back(
} (water_pressure_reservoir
- Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aquct_data_.d0))
* Base::alphai_[idx]);
}
// We take the average of the calculated equilibrium pressures. // We take the average of the calculated equilibrium pressures.
Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size(); Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / pw_aquifer.size();
return aquifer_pres_avg; return aquifer_pres_avg;
} }
}; // class AquiferCarterTracy }; // class AquiferCarterTracy
} // namespace Opm } // namespace Opm
#endif #endif

284
opm/simulators/aquifers/AquiferFetkovich.hpp
View File

@ -31,43 +31,44 @@ along with OPM. If not, see <http://www.gnu.org/licenses/>.
namespace Opm namespace Opm
{ {
template<typename TypeTag> template <typename TypeTag>
class AquiferFetkovich: public AquiferInterface<TypeTag> class AquiferFetkovich : public AquiferInterface<TypeTag>
{ {
public: public:
typedef AquiferInterface<TypeTag> Base; typedef AquiferInterface<TypeTag> Base;
using typename Base::Simulator;
using typename Base::ElementContext;
using typename Base::FluidSystem;
using typename Base::BlackoilIndices; using typename Base::BlackoilIndices;
using typename Base::RateVector; using typename Base::ElementContext;
using typename Base::IntensiveQuantities;
using typename Base::Eval; using typename Base::Eval;
using typename Base::Scalar;
using typename Base::FluidState; using typename Base::FluidState;
using typename Base::FluidSystem;
using typename Base::IntensiveQuantities;
using typename Base::RateVector;
using typename Base::Scalar;
using typename Base::Simulator;
using Base::waterCompIdx; using Base::waterCompIdx;
using Base::waterPhaseIdx; using Base::waterPhaseIdx;
AquiferFetkovich( const Aquancon::AquanconOutput& connection, AquiferFetkovich(const Aquancon::AquanconOutput& connection,
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(connection, 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& Qai : Base::Qai_) {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize(); Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize();
aquifer_pressure_ = aquiferPressure(); aquifer_pressure_ = aquiferPressure();
} }
} }
protected: protected:
// Aquifer Fetkovich Specific Variables // Aquifer Fetkovich Specific Variables
// TODO: using const reference here will cause segmentation fault, which is very strange // TODO: using const reference here will cause segmentation fault, which is very strange
const Aquifetp::AQUFETP_data aqufetp_data_; const Aquifetp::AQUFETP_data aqufetp_data_;
@ -75,173 +76,174 @@ namespace Opm
inline void initializeConnections(const Aquancon::AquanconOutput& connection) override inline void initializeConnections(const Aquancon::AquanconOutput& connection) override
{ {
const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid(); const auto& grid = eclState.getInputGrid();
Base::cell_idx_ = connection.global_index; Base::cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell(); auto globalCellIdx = ugrid.globalCell();
assert( Base::cell_idx_ == connection.global_index); assert(Base::cell_idx_ == connection.global_index);
assert( (Base::cell_idx_.size() == connection.influx_coeff.size()) ); assert((Base::cell_idx_.size() == connection.influx_coeff.size()));
assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); assert((connection.influx_coeff.size() == connection.influx_multiplier.size()));
assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); assert((connection.influx_multiplier.size() == connection.reservoir_face_dir.size()));
// 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); Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0);
Base::alphai_.resize(Base::cell_idx_.size(), 1.0); Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0); Base::faceArea_connected_.resize(Base::cell_idx_.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);
// Translate the C face tag into the enum used by opm-parser's TransMult class // Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection; Opm::FaceDir::DirEnum faceDirection;
// 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); Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) {
{ const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); Base::cellToConnectionIdx_[cell_index] = idx;
Base::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) {
{ // The index of the face in the compressed grid
// The index of the face in the compressed grid const int faceIdx = *cellFaceIter;
const int faceIdx = *cellFaceIter;
// the logically-Cartesian direction of the face // the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
switch(faceTag) switch (faceTag) {
{ case 0:
case 0: faceDirection = Opm::FaceDir::XMinus; faceDirection = Opm::FaceDir::XMinus;
break; break;
case 1: faceDirection = Opm::FaceDir::XPlus; case 1:
break; faceDirection = Opm::FaceDir::XPlus;
case 2: faceDirection = Opm::FaceDir::YMinus; break;
break; case 2:
case 3: faceDirection = Opm::FaceDir::YPlus; faceDirection = Opm::FaceDir::YMinus;
break; break;
case 4: faceDirection = Opm::FaceDir::ZMinus; case 3:
break; faceDirection = Opm::FaceDir::YPlus;
case 5: faceDirection = Opm::FaceDir::ZPlus; break;
break; case 4:
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer problem. Make sure faceTag is correctly defined"); faceDirection = Opm::FaceDir::ZMinus;
} break;
case 5:
faceDirection = Opm::FaceDir::ZPlus;
break;
default:
OPM_THROW(Opm::NumericalIssue,
"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
}
if (faceDirection == connection.reservoir_face_dir.at(idx)) if (faceDirection == connection.reservoir_face_dir.at(idx)) {
{ Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection);
Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); denom_face_areas += (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx));
denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) ); }
} }
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::cell_depth_.at(idx) = cellCenter[2];
} }
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::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 < Base::cell_idx_.size(); ++idx) {
{ Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)
Base::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.
: ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas; : (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas;
} }
} }
void assignRestartData(const data::AquiferData& xaq) override void assignRestartData(const data::AquiferData& xaq) override
{ {
if (xaq.type != data::AquiferType::Fetkovich) if (xaq.type != data::AquiferType::Fetkovich) {
{ throw std::invalid_argument {"Analytic aquifer data for unexpected aquifer type "
throw std::invalid_argument { "passed to Fetkovich aquifer"};
"Analytic aquifer data for unexpected aquifer type " }
"passed to Fetkovich aquifer"
};
}
this->aquifer_pressure_ = xaq.pressure; this->aquifer_pressure_ = xaq.pressure;
} }
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_ - Base::pressure_current_.at(idx)
+ Base::rhow_[idx] * Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0); + Base::rhow_[idx] * Base::gravity_() * (Base::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 = Base::W_flux_.value();
Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 ); Scalar pa_ = Base::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 ; Base::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() / Base::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; Base::Qai_.at(idx) = Base::alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef;
} }
inline void calculateAquiferCondition() override inline void calculateAquiferCondition() override
{ {
Base::rhow_.resize(Base::cell_idx_.size(),0.); Base::rhow_.resize(Base::cell_idx_.size(), 0.);
if (this->solution_set_from_restart_) { if (this->solution_set_from_restart_) {
return; return;
} }
if (!aqufetp_data_.p0) if (!aqufetp_data_.p0) {
{ Base::pa0_ = calculateReservoirEquilibrium();
Base::pa0_ = calculateReservoirEquilibrium(); } else {
} Base::pa0_ = *(aqufetp_data_.p0);
else }
{ aquifer_pressure_ = Base::pa0_;
Base::pa0_ = *(aqufetp_data_.p0);
}
aquifer_pressure_ = Base::pa0_ ;
} }
inline Scalar calculateReservoirEquilibrium() override inline Scalar calculateReservoirEquilibrium() override
{ {
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
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(Base::ebos_simulator_);
const auto& gridView = Base::ebos_simulator_.gridView(); const auto& gridView = Base::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 = Base::cellToConnectionIdx_[cellIdx];
if (idx < 0) if (idx < 0)
continue; continue;
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);
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); Base::rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back( (water_pressure_reservoir - Base::rhow_[idx].value()*Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0))*Base::alphai_[idx] ); pw_aquifer.push_back(
} (water_pressure_reservoir
- Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aqufetp_data_.d0))
* Base::alphai_[idx]);
}
// We take the average of the calculated equilibrium pressures. // We take the average of the calculated equilibrium pressures.
const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.); const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.);
const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha; const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha;
return aquifer_pres_avg; return aquifer_pres_avg;
} }
}; //Class AquiferFetkovich }; // Class AquiferFetkovich
} // namespace Opm } // namespace Opm
#endif #endif

223
opm/simulators/aquifers/AquiferInterface.hpp
View File

@ -22,28 +22,28 @@
#ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED #ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED #define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp> #include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifetp.hpp> #include <opm/parser/eclipse/EclipseState/Aquifetp.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/output/data/Aquifer.hpp> #include <opm/output/data/Aquifer.hpp>
#include <opm/material/common/MathToolbox.hpp> #include <opm/material/common/MathToolbox.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp> #include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp> #include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <vector>
#include <algorithm> #include <algorithm>
#include <unordered_map> #include <unordered_map>
#include <vector>
namespace Opm namespace Opm
{ {
template<typename TypeTag> template <typename TypeTag>
class AquiferInterface class AquiferInterface
{ {
public: public:
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@ -59,154 +59,165 @@ namespace Opm
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval; typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef Opm::BlackOilFluidState<Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, BlackoilIndices::numPhases> FluidState; typedef Opm::BlackOilFluidState<Eval,
FluidSystem,
enableTemperature,
enableEnergy,
BlackoilIndices::gasEnabled,
BlackoilIndices::numPhases>
FluidState;
static const auto waterCompIdx = FluidSystem::waterCompIdx; static const auto waterCompIdx = FluidSystem::waterCompIdx;
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
// Constructor // Constructor
AquiferInterface( const Aquancon::AquanconOutput& connection, AquiferInterface(const Aquancon::AquanconOutput& connection,
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) : connection_(connection)
, ebos_simulator_(ebosSimulator) , ebos_simulator_(ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed) , cartesian_to_compressed_(cartesian_to_compressed)
{} {
}
// Deconstructor // Deconstructor
virtual ~AquiferInterface() {} virtual ~AquiferInterface()
{
}
void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln) void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
{ {
auto xaqPos = std::find_if(aquiferSoln.begin(), aquiferSoln.end(), auto xaqPos
[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->connection_.aquiferID; });
});
if (xaqPos == aquiferSoln.end()) { if (xaqPos == aquiferSoln.end()) {
// No restart value applies to this aquifer. Nothing to do. // 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;
} }
void initialSolutionApplied() void initialSolutionApplied()
{ {
initQuantities(connection_); initQuantities(connection_);
} }
void beginTimeStep() void beginTimeStep()
{ {
ElementContext elemCtx(ebos_simulator_); ElementContext elemCtx(ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin<0>(); auto elemIt = ebos_simulator_.gridView().template begin<0>();
const auto& elemEndIt = ebos_simulator_.gridView().template end<0>(); const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt; const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryStencil(elem);
int cellIdx = elemCtx.globalSpaceIndex(0, 0); int cellIdx = elemCtx.globalSpaceIndex(0, 0);
int idx = cellToConnectionIdx_[cellIdx]; int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0) if (idx < 0)
continue; continue;
elemCtx.updateIntensiveQuantities(0); elemCtx.updateIntensiveQuantities(0);
const auto& iq = elemCtx.intensiveQuantities(0, 0); const auto& iq = elemCtx.intensiveQuantities(0, 0);
pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx));
} }
} }
template <class Context> template <class Context>
void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx)
{ {
unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int idx = cellToConnectionIdx_[cellIdx]; int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0) if (idx < 0)
return; return;
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const
// IntensiveQuantities of that particular cell_id // pointer to IntensiveQuantities of that particular cell_id
const IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx); const IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
// This is the pressure at td + dt // This is the pressure at td + dt
updateCellPressure(pressure_current_,idx,intQuants); updateCellPressure(pressure_current_, idx, intQuants);
updateCellDensity(idx,intQuants); updateCellDensity(idx, intQuants);
calculateInflowRate(idx, context.simulator()); calculateInflowRate(idx, context.simulator());
rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx]
Qai_[idx]/context.dofVolume(spaceIdx, timeIdx); += Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
} }
protected: protected:
inline Scalar gravity_() const inline Scalar gravity_() const
{ {
return ebos_simulator_.problem().gravity()[2]; return ebos_simulator_.problem().gravity()[2];
} }
inline void initQuantities(const Aquancon::AquanconOutput& connection) inline void initQuantities(const Aquancon::AquanconOutput& connection)
{ {
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0 // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
if (!this->solution_set_from_restart_) if (!this->solution_set_from_restart_) {
{ W_flux_ = 0.;
W_flux_ = 0.; }
}
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections function // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
initializeConnections(connection); // function
calculateAquiferCondition(); initializeConnections(connection);
calculateAquiferConstants(); calculateAquiferCondition();
calculateAquiferConstants();
pressure_previous_.resize(cell_idx_.size(), 0.); pressure_previous_.resize(cell_idx_.size(), 0.);
pressure_current_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.);
Qai_.resize(cell_idx_.size(), 0.0); Qai_.resize(cell_idx_.size(), 0.0);
} }
inline void updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants) inline void
updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{ {
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx); pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
} }
inline void updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants) inline void
updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{ {
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value();
} }
inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants)
{ {
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
rhow_.at(idx) = fs.density(waterPhaseIdx); rhow_.at(idx) = fs.density(waterPhaseIdx);
} }
template<class faceCellType, class ugridType> template <class faceCellType, class ugridType>
inline double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid, inline double getFaceArea(const faceCellType& faceCells,
const int faceIdx, const int idx, const ugridType& ugrid,
const int faceIdx,
const int idx,
const Aquancon::AquanconOutput& connection) const const Aquancon::AquanconOutput& connection) const
{ {
// Check now if the face is outside of the reservoir, or if it adjoins an inactive cell // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell
// Do not make the connection if the product of the two cellIdx > 0. This is because the // Do not make the connection if the product of the two cellIdx > 0. This is because the
// face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell
double faceArea = 0.; // adjoining)
const auto cellNeighbour0 = faceCells(faceIdx,0); double faceArea = 0.;
const auto cellNeighbour1 = faceCells(faceIdx,1); const auto cellNeighbour0 = faceCells(faceIdx, 0);
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); const auto cellNeighbour1 = faceCells(faceIdx, 1);
const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))? const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
defaultFaceArea : const auto calculatedFaceArea
*(connection.influx_coeff.at(idx)); = (!connection.influx_coeff.at(idx)) ? defaultFaceArea : *(connection.influx_coeff.at(idx));
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.;
} } else if (cellNeighbour0 == 0) {
else if (cellNeighbour0 == 0){ faceArea = (cellNeighbour1 < 0) ? faceArea : 0.;
faceArea = (cellNeighbour1 < 0)? faceArea : 0.; }
} return faceArea;
return faceArea;
} }
virtual void endTimeStep() = 0; virtual void endTimeStep() = 0;
@ -232,9 +243,9 @@ namespace Opm
Eval W_flux_; Eval W_flux_;
bool solution_set_from_restart_{false}; bool solution_set_from_restart_ {false};
virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; virtual void initializeConnections(const Aquancon::AquanconOutput& connection) = 0;
virtual void assignRestartData(const data::AquiferData& xaq) = 0; virtual void assignRestartData(const data::AquiferData& xaq) = 0;
@ -244,8 +255,8 @@ namespace Opm
virtual void calculateAquiferConstants() = 0; virtual void calculateAquiferConstants() = 0;
virtual Scalar calculateReservoirEquilibrium() =0; virtual Scalar calculateReservoirEquilibrium() = 0;
// 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
}; };
} // namespace Opm } // namespace Opm
#endif #endif

87
opm/simulators/aquifers/BlackoilAquiferModel.hpp
View File

@ -26,9 +26,9 @@
#include <ebos/eclbaseaquifermodel.hh> #include <ebos/eclbaseaquifermodel.hh>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp> #include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifetp.hpp> #include <opm/parser/eclipse/EclipseState/Aquifetp.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/output/data/Aquifer.hpp> #include <opm/output/data/Aquifer.hpp>
@ -39,62 +39,59 @@
#include <vector> #include <vector>
namespace Opm { namespace Opm
{
/// Class for handling the blackoil well model. /// Class for handling the blackoil well model.
template<typename TypeTag> template <typename TypeTag>
class BlackoilAquiferModel class BlackoilAquiferModel
{ {
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
public: public:
explicit BlackoilAquiferModel(Simulator& simulator); explicit BlackoilAquiferModel(Simulator& simulator);
void initialSolutionApplied(); void initialSolutionApplied();
void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln); void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln);
void beginEpisode(); void beginEpisode();
void beginTimeStep(); void beginTimeStep();
void beginIteration(); void beginIteration();
// add the water rate due to aquifers to the source term. // add the water rate due to aquifers to the source term.
template <class Context> template <class Context>
void addToSource(RateVector& rates, void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const;
const Context& context, void endIteration();
unsigned spaceIdx, void endTimeStep();
unsigned timeIdx) const; void endEpisode();
void endIteration();
void endTimeStep();
void endEpisode();
template <class Restarter> template <class Restarter>
void serialize(Restarter& res); void serialize(Restarter& res);
template <class Restarter> template <class Restarter>
void deserialize(Restarter& res); void deserialize(Restarter& res);
protected: protected:
// --------- Types --------- // --------- Types ---------
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef AquiferCarterTracy<TypeTag> AquiferCarterTracy_object; typedef AquiferCarterTracy<TypeTag> AquiferCarterTracy_object;
typedef AquiferFetkovich<TypeTag> AquiferFetkovich_object; typedef AquiferFetkovich<TypeTag> AquiferFetkovich_object;
Simulator& simulator_; Simulator& simulator_;
std::unordered_map<int, int> cartesian_to_compressed_; std::unordered_map<int, int> cartesian_to_compressed_;
mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy; mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy;
mutable std::vector<AquiferFetkovich_object> aquifers_Fetkovich; mutable std::vector<AquiferFetkovich_object> aquifers_Fetkovich;
// This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
void init(); void init();
bool aquiferActive() const; bool aquiferActive() const;
bool aquiferCarterTracyActive() const; bool aquiferCarterTracyActive() const;
bool aquiferFetkovichActive() const; bool aquiferFetkovichActive() const;
};
};
} // namespace Opm } // namespace Opm

349
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
View File

@ -1,225 +1,208 @@
#include <opm/grid/utility/cartesianToCompressed.hpp> #include <opm/grid/utility/cartesianToCompressed.hpp>
namespace Opm { namespace Opm
{
template<typename TypeTag> template <typename TypeTag>
BlackoilAquiferModel<TypeTag>:: BlackoilAquiferModel<TypeTag>::BlackoilAquiferModel(Simulator& simulator)
BlackoilAquiferModel(Simulator& simulator) : simulator_(simulator)
: simulator_(simulator) {
{
init(); init();
} }
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::initialSolutionApplied() BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
{ {
if(aquiferCarterTracyActive()) if (aquiferCarterTracyActive()) {
{ for (auto& aquifer : aquifers_CarterTracy) {
for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) aquifer.initialSolutionApplied();
{ }
aquifer->initialSolutionApplied();
}
} }
if(aquiferFetkovichActive()) if (aquiferFetkovichActive()) {
{ for (auto& aquifer : aquifers_Fetkovich) {
for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) aquifer.initialSolutionApplied();
{ }
aquifer->initialSolutionApplied();
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferData>& aquiferSoln) BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
{ {
if(aquiferCarterTracyActive()) if (aquiferCarterTracyActive()) {
{ for (auto& aquifer : aquifers_CarterTracy) {
for (auto& aquifer : aquifers_CarterTracy) aquifer.initFromRestart(aquiferSoln);
{ }
aquifer.initFromRestart(aquiferSoln);
}
} }
if(aquiferFetkovichActive()) if (aquiferFetkovichActive()) {
{ for (auto& aquifer : aquifers_Fetkovich) {
for (auto& aquifer : aquifers_Fetkovich) aquifer.initFromRestart(aquiferSoln);
{ }
aquifer.initFromRestart(aquiferSoln);
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::beginEpisode() BlackoilAquiferModel<TypeTag>::beginEpisode()
{ } {
}
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::beginIteration() BlackoilAquiferModel<TypeTag>::beginIteration()
{ } {
}
template<typename TypeTag> template <typename TypeTag>
void BlackoilAquiferModel<TypeTag>:: beginTimeStep() void
{ BlackoilAquiferModel<TypeTag>::beginTimeStep()
if(aquiferCarterTracyActive()) {
{ if (aquiferCarterTracyActive()) {
for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) for (auto& aquifer : aquifers_CarterTracy) {
{ aquifer.beginTimeStep();
aquifer->beginTimeStep(); }
}
} }
if(aquiferFetkovichActive()) if (aquiferFetkovichActive()) {
{ for (auto& aquifer : aquifers_Fetkovich) {
for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) aquifer.beginTimeStep();
{ }
aquifer->beginTimeStep();
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
template<class Context> template <class Context>
void BlackoilAquiferModel<TypeTag>:: addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const void
{ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
if(aquiferCarterTracyActive()) const Context& context,
{ unsigned spaceIdx,
for (auto& aquifer : aquifers_CarterTracy) unsigned timeIdx) const
{ {
aquifer.addToSource(rates, context, spaceIdx, timeIdx); if (aquiferCarterTracyActive()) {
} for (auto& aquifer : aquifers_CarterTracy) {
aquifer.addToSource(rates, context, spaceIdx, timeIdx);
}
} }
if(aquiferFetkovichActive()) if (aquiferFetkovichActive()) {
{ for (auto& aquifer : aquifers_Fetkovich) {
for (auto& aquifer : aquifers_Fetkovich) aquifer.addToSource(rates, context, spaceIdx, timeIdx);
{ }
aquifer.addToSource(rates, context, spaceIdx, timeIdx);
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::endIteration() BlackoilAquiferModel<TypeTag>::endIteration()
{ } {
}
template<typename TypeTag> template <typename TypeTag>
void BlackoilAquiferModel<TypeTag>:: endTimeStep() void
{ BlackoilAquiferModel<TypeTag>::endTimeStep()
if(aquiferCarterTracyActive()) {
{ if (aquiferCarterTracyActive()) {
for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) for (auto& aquifer : aquifers_CarterTracy) {
{ aquifer.endTimeStep();
aquifer->endTimeStep(); }
}
} }
if(aquiferFetkovichActive()) if (aquiferFetkovichActive()) {
{ for (auto& aquifer : aquifers_Fetkovich) {
for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) aquifer.endTimeStep();
{ }
aquifer->endTimeStep();
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>::endEpisode() BlackoilAquiferModel<TypeTag>::endEpisode()
{ } {
}
template <typename TypeTag> template <typename TypeTag>
template <class Restarter> template <class Restarter>
void void
BlackoilAquiferModel<TypeTag>::serialize(Restarter& /* res */) BlackoilAquiferModel<TypeTag>::serialize(Restarter& /* res */)
{ {
// TODO (?) // TODO (?)
throw std::logic_error("BlackoilAquiferModel::serialize() is not yet implemented"); throw std::logic_error("BlackoilAquiferModel::serialize() is not yet implemented");
} }
template<typename TypeTag> template <typename TypeTag>
template <class Restarter> template <class Restarter>
void void
BlackoilAquiferModel<TypeTag>::deserialize(Restarter& /* res */) BlackoilAquiferModel<TypeTag>::deserialize(Restarter& /* res */)
{ {
// TODO (?) // TODO (?)
throw std::logic_error("BlackoilAquiferModel::deserialize() is not yet implemented"); throw std::logic_error("BlackoilAquiferModel::deserialize() is not yet implemented");
} }
// Initialize the aquifers in the deck // Initialize the aquifers in the deck
template<typename TypeTag> template <typename TypeTag>
void void
BlackoilAquiferModel<TypeTag>:: init() BlackoilAquiferModel<TypeTag>::init()
{ {
const auto& deck = this->simulator_.vanguard().deck(); const auto& deck = this->simulator_.vanguard().deck();
if (deck.hasKeyword("AQUCT")) { if (deck.hasKeyword("AQUCT")) {
//updateConnectionIntensiveQuantities(); // updateConnectionIntensiveQuantities();
const auto& eclState = this->simulator_.vanguard().eclState(); const auto& eclState = this->simulator_.vanguard().eclState();
// Get all the carter tracy aquifer properties data and put it in aquifers vector // Get all the carter tracy aquifer properties data and put it in aquifers vector
const AquiferCT aquiferct = AquiferCT(eclState,deck); const AquiferCT aquiferct = AquiferCT(eclState, deck);
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.getAquifers(); std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.getAquifers();
std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput(); std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
assert( aquifersData.size() == aquifer_connection.size() ); 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, cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
Opm::UgGridHelpers::globalCell(ugrid));
for (size_t i = 0; i < aquifersData.size(); ++i) for (size_t i = 0; i < aquifersData.size(); ++i) {
{ aquifers_CarterTracy.push_back(AquiferCarterTracy<TypeTag>(
aquifers_CarterTracy.push_back( aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_, aquifersData.at(i)));
AquiferCarterTracy<TypeTag> (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) }
);
}
} }
if(deck.hasKeyword("AQUFETP")) if (deck.hasKeyword("AQUFETP")) {
{ // updateConnectionIntensiveQuantities();
//updateConnectionIntensiveQuantities(); const auto& eclState = this->simulator_.vanguard().eclState();
const auto& eclState = this->simulator_.vanguard().eclState();
// Get all the carter tracy aquifer properties data and put it in aquifers vector // Get all the carter tracy aquifer properties data and put it in aquifers vector
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.getAquifers(); std::vector<Aquifetp::AQUFETP_data> aquifersData = aquifetp.getAquifers();
std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput(); std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
assert( aquifersData.size() == aquifer_connection.size() ); 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, cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid));
Opm::UgGridHelpers::globalCell(ugrid));
for (size_t i = 0; i < aquifersData.size(); ++i) for (size_t i = 0; i < aquifersData.size(); ++i) {
{ aquifers_Fetkovich.push_back(AquiferFetkovich<TypeTag>(
aquifers_Fetkovich.push_back( aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_, aquifersData.at(i)));
AquiferFetkovich<TypeTag> (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) }
);
}
} }
} }
template<typename TypeTag> template <typename TypeTag>
bool bool
BlackoilAquiferModel<TypeTag>:: aquiferActive() const BlackoilAquiferModel<TypeTag>::aquiferActive() const
{ {
return (aquiferCarterTracyActive() || aquiferFetkovichActive()); return (aquiferCarterTracyActive() || aquiferFetkovichActive());
} }
template<typename TypeTag> template <typename TypeTag>
bool bool
BlackoilAquiferModel<TypeTag>:: aquiferCarterTracyActive() const BlackoilAquiferModel<TypeTag>::aquiferCarterTracyActive() const
{ {
return !aquifers_CarterTracy.empty(); return !aquifers_CarterTracy.empty();
} }
template<typename TypeTag> template <typename TypeTag>
bool bool
BlackoilAquiferModel<TypeTag>:: aquiferFetkovichActive() const BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
{ {
return !aquifers_Fetkovich.empty(); return !aquifers_Fetkovich.empty();
} }
} // namespace Opm } // namespace Opm