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 #1695 from WesselDeZeeuw/BaseClassAquifer

Browse Source 
Base class aquifer
This commit is contained in:
Kai Bao 2019-02-18 10:44:51 +01:00 committed by GitHub
commit 2d933f2303
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 572 additions and 655 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -104,7 +104,9 @@ list (APPEND TEST_DATA_FILES
# originally generated with the command: # originally generated with the command:
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
list (APPEND PUBLIC_HEADER_FILES list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AquiferInterface.hpp
opm/autodiff/AquiferCarterTracy.hpp opm/autodiff/AquiferCarterTracy.hpp
opm/autodiff/AquiferFetkovich.hpp
opm/autodiff/BlackoilAmg.hpp opm/autodiff/BlackoilAmg.hpp
opm/autodiff/BlackoilDetails.hpp opm/autodiff/BlackoilDetails.hpp
opm/autodiff/BlackoilModelParametersEbos.hpp opm/autodiff/BlackoilModelParametersEbos.hpp

368
opm/autodiff/AquiferCarterTracy.hpp
View File

@ -21,268 +21,68 @@
#ifndef OPM_AQUIFERCT_HEADER_INCLUDED #ifndef OPM_AQUIFERCT_HEADER_INCLUDED
#define OPM_AQUIFERCT_HEADER_INCLUDED #define OPM_AQUIFERCT_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp> #include <opm/autodiff/AquiferInterface.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp>
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <vector>
#include <algorithm>
#include <unordered_map>
namespace Opm namespace Opm
{ {
template<typename TypeTag> template<typename TypeTag>
class AquiferCarterTracy class AquiferCarterTracy: public AquiferInterface<TypeTag>
{ {
public: public:
typedef AquiferInterface<TypeTag> Base;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; using typename Base::Simulator;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; using typename Base::ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; using typename Base::FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; using typename Base::BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; using typename Base::RateVector;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; using typename Base::IntensiveQuantities;
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; using typename Base::Eval;
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; using typename Base::Scalar;
using typename Base::FluidState;
static const int numEq = BlackoilIndices::numEq; using Base::waterCompIdx;
typedef double Scalar; using Base::waterPhaseIdx;
AquiferCarterTracy( const Aquancon::AquanconOutput& connection,
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef Opm::BlackOilFluidState<Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, BlackoilIndices::numPhases> FluidState;
static const auto waterCompIdx = FluidSystem::waterCompIdx;
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data,
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,
: ebos_simulator_ (ebosSimulator) const AquiferCT::AQUCT_data& aquct_data)
, aquct_data_ (aquct_data) : Base(connection, cartesian_to_compressed, ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed) , aquct_data_(aquct_data)
, connection_(connection)
{} {}
void initialSolutionApplied()
{
initQuantities(connection_);
}
void beginTimeStep()
{
ElementContext elemCtx(ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin<0>();
const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
int cellIdx = elemCtx.globalSpaceIndex(0, 0);
int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updateIntensiveQuantities(0);
const auto& iq = elemCtx.intensiveQuantities(0, 0);
pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx));
}
}
template <class Context>
void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx)
{
unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
return;
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to
// IntensiveQuantities of that particular cell_id
const IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
// This is the pressure at td + dt
updateCellPressure(pressure_current_,idx,intQuants);
updateCellDensity(idx,intQuants);
calculateInflowRate(idx, context.simulator());
rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] +=
Qai_[idx]/context.dofVolume(spaceIdx, timeIdx);
}
void endTimeStep() void endTimeStep()
{ {
for (const auto& Qai: Qai_) { for (const auto& Qai: Base::Qai_) {
W_flux_ += Qai*ebos_simulator_.timeStepSize(); Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize();
} }
} }
private: protected:
const Simulator& ebos_simulator_;
// Grid variables
std::vector<size_t> cell_idx_;
std::vector<Scalar> faceArea_connected_;
// Quantities at each grid id
std::vector<Scalar> cell_depth_;
std::vector<Scalar> pressure_previous_;
std::vector<Eval> pressure_current_;
std::vector<Eval> Qai_;
std::vector<Eval> rhow_;
std::vector<Scalar> alphai_;
// Variables constants // Variables constants
const AquiferCT::AQUCT_data aquct_data_; const AquiferCT::AQUCT_data aquct_data_;
Scalar mu_w_; //water viscosity
Scalar beta_; // Influx constant Scalar beta_; // Influx constant
Scalar Tc_; // Time constant
Scalar pa0_; // initial aquifer pressure
Eval W_flux_;
const std::unordered_map<int, int>& cartesian_to_compressed_;
Scalar gravity_() const
{ return ebos_simulator_.problem().gravity()[2]; }
inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
{
// We use the opm-common numeric linear interpolator
pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td);
pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td);
}
inline void initQuantities(const Aquancon::AquanconOutput& connection)
{
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0
W_flux_ = 0.;
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections function
initializeConnections(connection);
calculateAquiferCondition();
calculateAquiferConstants();
pressure_previous_.resize(cell_idx_.size(), 0.);
pressure_current_.resize(cell_idx_.size(), 0.);
Qai_.resize(cell_idx_.size(), 0.0);
}
inline void updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
}
inline void updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value();
}
inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
rhow_.at(idx) = fs.density(waterPhaseIdx);
}
inline Scalar dpai(int idx)
{
Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_()*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx);
return dp;
}
// 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)
{
const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Tc_;
const Scalar td = simulator.time() / Tc_;
Scalar PItdprime = 0.;
Scalar PItd = 0.;
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime );
b = beta_ / (Tc_ * ( PItd - td*PItdprime));
}
// This function implements Eq 5.7 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator)
{
Scalar a, b;
calculateEqnConstants(a,b,idx,simulator);
Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) );
}
inline void calculateAquiferConstants()
{
// 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
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 );
}
template<class faceCellType, class ugridType>
inline const double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid,
const int faceIdx, const int idx,
const Aquancon::AquanconOutput& connection) const
{
// 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
// face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining)
double faceArea = 0.;
const auto cellNeighbour0 = faceCells(faceIdx,0);
const auto cellNeighbour1 = faceCells(faceIdx,1);
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))?
defaultFaceArea :
*(connection.influx_coeff.at(idx));
faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea;
if (cellNeighbour1 == 0){
faceArea = (cellNeighbour0 < 0)? faceArea : 0.;
}
else if (cellNeighbour0 == 0){
faceArea = (cellNeighbour1 < 0)? faceArea : 0.;
}
return faceArea;
}
// 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(const Aquancon::AquanconOutput& connection) inline void initializeConnections(const Aquancon::AquanconOutput& connection)
{ {
const auto& eclState = ebos_simulator_.vanguard().eclState(); const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
const auto& ugrid = ebos_simulator_.vanguard().grid(); const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid(); const auto& grid = eclState.getInputGrid();
cell_idx_ = connection.global_index; Base::cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell(); auto globalCellIdx = ugrid.globalCell();
assert( cell_idx_ == connection.global_index); assert( Base::cell_idx_ == connection.global_index);
assert( (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
cell_depth_.resize(cell_idx_.size(), aquct_data_.d0); Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0);
alphai_.resize(cell_idx_.size(), 1.0); Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
faceArea_connected_.resize(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);
@ -292,11 +92,11 @@ namespace Opm
// 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.;
cellToConnectionIdx_.resize(ebos_simulator_.gridView().size(/*codim=*/0), -1); Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < cell_idx_.size(); ++idx) for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx)
{ {
const int cell_index = cartesian_to_compressed_.at(cell_idx_[idx]); const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
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)
@ -326,48 +126,91 @@ namespace Opm
if (faceDirection == connection.reservoir_face_dir.at(idx)) if (faceDirection == connection.reservoir_face_dir.at(idx))
{ {
faceArea_connected_.at(idx) = 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) * faceArea_connected_.at(idx) ); denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) );
} }
} }
auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
cell_depth_.at(idx) = cellCenter[2]; 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 < cell_idx_.size(); ++idx)
{ {
alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
0. 0.
: ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas; : ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas;
} }
} }
inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
{
// We use the opm-common numeric linear interpolator
pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td);
pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td);
}
inline Scalar dpai(int idx)
{
Scalar dp = Base::pa0_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aquct_data_.d0) - Base::pressure_previous_.at(idx);
return dp;
}
// 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)
{
const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_;
const Scalar td = simulator.time() / Base::Tc_;
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));
}
// This function implements Eq 5.7 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator)
{
Scalar a, b;
calculateEqnConstants(a,b,idx,simulator);
Base::Qai_.at(idx) = Base::alphai_.at(idx)*( a - b * ( Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx) ) );
}
inline void calculateAquiferConstants()
{
// 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_ = Base::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() inline void calculateAquiferCondition()
{ {
int pvttableIdx = aquct_data_.pvttableID - 1; int pvttableIdx = aquct_data_.pvttableID - 1;
Base::rhow_.resize(Base::cell_idx_.size(),0.);
rhow_.resize(cell_idx_.size(),0.);
if (!aquct_data_.p0) if (!aquct_data_.p0)
{ {
pa0_ = calculateReservoirEquilibrium(); Base::pa0_ = calculateReservoirEquilibrium();
} }
else else
{ {
pa0_ = *(aquct_data_.p0); Base::pa0_ = *(aquct_data_.p0);
} }
// 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(ebos_simulator_); ElementContext elemCtx(Base::ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin</*codim=*/0>(); auto elemIt = Base::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);
// Initialize a FluidState object first // Initialize a FluidState object first
FluidState fs_aquifer; FluidState fs_aquifer;
// We use the temperature of the first cell connected to the aquifer // We use the temperature of the first cell connected to the aquifer
@ -375,11 +218,9 @@ namespace Opm
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 = pa0_; pa0_mean = Base::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);
Base::mu_w_ = mu_w_aquifer.value();
mu_w_ = mu_w_aquifer.value();
} }
@ -390,8 +231,8 @@ namespace Opm
std::vector<Scalar> pw_aquifer; std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir; Scalar water_pressure_reservoir;
ElementContext elemCtx(ebos_simulator_); ElementContext elemCtx(Base::ebos_simulator_);
const auto& gridView = 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) {
@ -399,29 +240,24 @@ namespace Opm
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 = 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();
rhow_[idx] = fs.density(waterPhaseIdx); Base::rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back( (water_pressure_reservoir - rhow_[idx].value()*gravity_()*(cell_depth_[idx] - aquct_data_.d0))*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;
} }
const Aquancon::AquanconOutput connection_;
std::vector<int> cellToConnectionIdx_;
}; // class AquiferCarterTracy }; // class AquiferCarterTracy
} // namespace Opm } // namespace Opm
#endif #endif

627
opm/autodiff/AquiferFetkovich.hpp
View File

@ -1,387 +1,240 @@
/* /*
Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface
Copyright 2017 Statoil ASA. Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or the Free Software Foundation, either version 3 of the License, or
(at your option) any later version. (at your option) any later version.
OPM is distributed in the hope that it will be useful, OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef OPM_AQUIFETP_HEADER_INCLUDED #ifndef OPM_AQUIFETP_HEADER_INCLUDED
#define OPM_AQUIFETP_HEADER_INCLUDED #define OPM_AQUIFETP_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Aquifetp.hpp> #include <opm/autodiff/AquiferInterface.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp> namespace Opm
#include <opm/common/utility/numeric/linearInterpolation.hpp> {
#include <opm/material/common/MathToolbox.hpp> template<typename TypeTag>
#include <opm/material/densead/Math.hpp> class AquiferFetkovich: public AquiferInterface<TypeTag>
#include <opm/material/densead/Evaluation.hpp> {
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
public:
#include <vector> typedef AquiferInterface<TypeTag> Base;
#include <algorithm>
#include <unordered_map> using typename Base::Simulator;
using typename Base::ElementContext;
namespace Opm using typename Base::FluidSystem;
{ using typename Base::BlackoilIndices;
using typename Base::RateVector;
template<typename TypeTag> using typename Base::IntensiveQuantities;
class AquiferFetkovich using typename Base::Eval;
{ using typename Base::Scalar;
using typename Base::FluidState;
public:
using Base::waterCompIdx;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; using Base::waterPhaseIdx;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; AquiferFetkovich( const Aquancon::AquanconOutput& connection,
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; const std::unordered_map<int, int>& cartesian_to_compressed,
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; const Simulator& ebosSimulator,
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; const Aquifetp::AQUFETP_data& aqufetp_data)
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; : Base(connection, cartesian_to_compressed, ebosSimulator)
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; , aqufetp_data_(aqufetp_data)
{}
static const int numEq = BlackoilIndices::numEq;
typedef double Scalar; void endTimeStep()
{
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval; for (const auto& Qai: Base::Qai_) {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize();
typedef Opm::BlackOilFluidState<Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, BlackoilIndices::numPhases> FluidState; aquifer_pressure_ = aquiferPressure();
}
static const auto waterCompIdx = FluidSystem::waterCompIdx; }
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
protected:
AquiferFetkovich( const Aquifetp::AQUFETP_data& aqufetp_data, // Aquifer Fetkovich Specific Variables
const Aquancon::AquanconOutput& connection, const Aquifetp::AQUFETP_data aqufetp_data_;
const std::unordered_map<int, int>& cartesian_to_compressed, Scalar aquifer_pressure_; // aquifer
const Simulator& ebosSimulator)
: ebos_simulator_ (ebosSimulator) inline void initializeConnections(const Aquancon::AquanconOutput& connection)
, cartesian_to_compressed_(cartesian_to_compressed) {
, aqufetp_data_ (aqufetp_data) const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
, connection_ (connection) const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
{} const auto& grid = eclState.getInputGrid();
void initialSolutionApplied() Base::cell_idx_ = connection.global_index;
{ auto globalCellIdx = ugrid.globalCell();
initQuantities(connection_);
} assert( Base::cell_idx_ == connection.global_index);
assert( (Base::cell_idx_.size() == connection.influx_coeff.size()) );
void beginTimeStep() assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) );
{ assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) );
ElementContext elemCtx(ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin<0>(); // We hack the cell depth values for now. We can actually get it from elementcontext pos
const auto& elemEndIt = ebos_simulator_.gridView().template end<0>(); Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0);
for (; elemIt != elemEndIt; ++elemIt) { Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
const auto& elem = *elemIt; Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0);
elemCtx.updatePrimaryStencil(elem); auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
int cellIdx = elemCtx.globalSpaceIndex(0, 0);
int idx = cellToConnectionIdx_[cellIdx]; // Translate the C face tag into the enum used by opm-parser's TransMult class
if (idx < 0) Opm::FaceDir::DirEnum faceDirection;
continue;
// denom_face_areas is the sum of the areas connected to an aquifer
elemCtx.updateIntensiveQuantities(0); Scalar denom_face_areas = 0.;
const auto& iq = elemCtx.intensiveQuantities(0, 0); Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); 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;
template <class Context>
void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const auto cellFacesRange = cell2Faces[cell_index];
{ for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); {
// The index of the face in the compressed grid
int idx = cellToConnectionIdx_[cellIdx]; const int faceIdx = *cellFaceIter;
if (idx < 0)
return; // the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to
// IntensiveQuantities of that particular cell_id switch(faceTag)
const IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx); {
// This is the pressure at td + dt case 0: faceDirection = Opm::FaceDir::XMinus;
updateCellPressure(pressure_current_,idx,intQuants); break;
updateCellDensity(idx,intQuants); case 1: faceDirection = Opm::FaceDir::XPlus;
calculateInflowRate(idx, context.simulator()); break;
rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += case 2: faceDirection = Opm::FaceDir::YMinus;
Qai_[idx]/context.dofVolume(spaceIdx, timeIdx); break;
} case 3: faceDirection = Opm::FaceDir::YPlus;
break;
void endTimeStep() case 4: faceDirection = Opm::FaceDir::ZMinus;
{ break;
for (const auto& Qai: Qai_) { case 5: faceDirection = Opm::FaceDir::ZPlus;
W_flux_ += Qai*ebos_simulator_.timeStepSize(); break;
aquifer_pressure_ = aquiferPressure(); default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
} }
}
private: if (faceDirection == connection.reservoir_face_dir.at(idx))
const Simulator& ebos_simulator_; {
const std::unordered_map<int, int>& cartesian_to_compressed_; 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) );
// Grid variables }
std::vector<size_t> cell_idx_; }
std::vector<Scalar> faceArea_connected_; auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::cell_depth_.at(idx) = cellCenter[2];
// Quantities at each grid id }
std::vector<Scalar> cell_depth_;
std::vector<Scalar> pressure_previous_; const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
std::vector<Eval> pressure_current_; for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx)
std::vector<Eval> Qai_; {
std::vector<Eval> rhow_; Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
std::vector<Scalar> alphai_; 0.
std::vector<int> cellToConnectionIdx_; : ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas;
}
// Variables constants }
const Aquifetp::AQUFETP_data aqufetp_data_;
const Aquancon::AquanconOutput connection_; inline Scalar dpai(int idx)
{
Scalar mu_w_; //water viscosity Scalar dp = aquifer_pressure_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0) - Base::pressure_current_.at(idx).value() ;
Scalar Tc_; // Time Constant return dp;
Scalar pa0_; // initial aquifer pressure }
Scalar aquifer_pressure_; // aquifer pressure
// This function implements Eq 5.12 of the EclipseTechnicalDescription
Eval W_flux_; inline Scalar aquiferPressure()
{
Scalar gravity_() const Scalar Flux = Base::W_flux_.value();
{ return ebos_simulator_.problem().gravity()[2]; } Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 );
return pa_;
inline void initQuantities(const Aquancon::AquanconOutput& connection) }
{
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0 inline void calculateAquiferConstants()
W_flux_ = 0.; {
Base::Tc_ = ( aqufetp_data_.C_t * aqufetp_data_.V0 ) / aqufetp_data_.J ;
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections function }
initializeConnections(connection); // This function implements Eq 5.14 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator)
calculateAquiferCondition(); {
Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ;
pressure_previous_.resize(cell_idx_.size(), 0.); Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_;
pressure_current_.resize(cell_idx_.size(), 0.); Base::Qai_.at(idx) = Base::alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_;
Qai_.resize(cell_idx_.size(), 0.0); }
}
inline void calculateAquiferCondition()
inline void updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants) {
{ int pvttableIdx = aqufetp_data_.pvttableID - 1;
const auto& fs = intQuants.fluidState(); Base::rhow_.resize(Base::cell_idx_.size(),0.);
pressure_water.at(idx) = fs.pressure(waterPhaseIdx); if (!aqufetp_data_.p0)
} {
Base::pa0_ = calculateReservoirEquilibrium();
inline void updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants) }
{ else
const auto& fs = intQuants.fluidState(); {
pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); Base::pa0_ = *(aqufetp_data_.p0);
} }
aquifer_pressure_ = Base::pa0_ ;
inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants)
{ // use the thermodynamic state of the first active cell as a
const auto& fs = intQuants.fluidState(); // reference. there might be better ways to do this...
rhow_.at(idx) = fs.density(waterPhaseIdx); ElementContext elemCtx(Base::ebos_simulator_);
} auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>();
elemCtx.updatePrimaryStencil(*elemIt);
inline Scalar dpai(int idx) elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
{ const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
Scalar dp = aquifer_pressure_ + rhow_.at(idx).value()*gravity_()*(cell_depth_.at(idx) - aqufetp_data_.d0) - pressure_current_.at(idx).value() ; // Initialize a FluidState object first
return dp; FluidState fs_aquifer;
} // We use the temperature of the first cell connected to the aquifer
// This function implements Eq 5.12 of the EclipseTechnicalDescription // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
inline Scalar aquiferPressure() fs_aquifer.assign( iq0.fluidState() );
{ Eval temperature_aq, pa0_mean;
Scalar Flux = W_flux_.value(); temperature_aq = fs_aquifer.temperature(0);
Scalar pa_ = pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 ); pa0_mean = Base::pa0_;
return pa_; Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
} Base::mu_w_ = mu_w_aquifer.value();
// This function implements Eq 5.14 of the EclipseTechnicalDescription }
inline void calculateInflowRate(int idx, const Simulator& simulator)
{ inline Scalar calculateReservoirEquilibrium()
Tc_ = ( aqufetp_data_.C_t * aqufetp_data_.V0 ) / aqufetp_data_.J ; {
Scalar td_Tc_ = simulator.timeStepSize() / Tc_ ; // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_; std::vector<Scalar> pw_aquifer;
Qai_.at(idx) = alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_; Scalar water_pressure_reservoir;
}
ElementContext elemCtx(Base::ebos_simulator_);
template<class faceCellType, class ugridType> const auto& gridView = Base::ebos_simulator_.gridView();
inline const double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid, auto elemIt = gridView.template begin</*codim=*/0>();
const int faceIdx, const int idx, const auto& elemEndIt = gridView.template end</*codim=*/0>();
const Aquancon::AquanconOutput& connection) const for (; elemIt != elemEndIt; ++elemIt) {
{ const auto& elem = *elemIt;
// Check now if the face is outside of the reservoir, or if it adjoins an inactive cell elemCtx.updatePrimaryStencil(elem);
// Do not make the connection if the product of the two cellIdx > 0. This is because the size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
// face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) int idx = Base::cellToConnectionIdx_[cellIdx];
double faceArea = 0.; if (idx < 0)
const auto cellNeighbour0 = faceCells(faceIdx,0); continue;
const auto cellNeighbour1 = faceCells(faceIdx,1);
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))? const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
defaultFaceArea : const auto& fs = iq0.fluidState();
*(connection.influx_coeff.at(idx));
faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea; water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
if (cellNeighbour1 == 0){ Base::rhow_[idx] = fs.density(waterPhaseIdx);
faceArea = (cellNeighbour0 < 0)? faceArea : 0.; pw_aquifer.push_back( (water_pressure_reservoir - Base::rhow_[idx].value()*Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0))*Base::alphai_[idx] );
} }
else if (cellNeighbour0 == 0){
faceArea = (cellNeighbour1 < 0)? faceArea : 0.; // 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();
return faceArea; return aquifer_pres_avg;
} }
}; //Class AquiferFetkovich
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer } // namespace Opm
inline void initializeConnections(const Aquancon::AquanconOutput& connection) #endif
{
const auto& eclState = ebos_simulator_.vanguard().eclState();
const auto& ugrid = ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell();
assert( cell_idx_ == connection.global_index);
assert( (cell_idx_.size() == connection.influx_coeff.size()) );
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
cell_depth_.resize(cell_idx_.size(), aqufetp_data_.d0);
alphai_.resize(cell_idx_.size(), 1.0);
faceArea_connected_.resize(cell_idx_.size(),0.0);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas = 0.;
cellToConnectionIdx_.resize(ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
{
const int cell_index = cartesian_to_compressed_.at(cell_idx_[idx]);
cellToConnectionIdx_[cell_index] = idx;
const auto cellFacesRange = cell2Faces[cell_index];
for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
{
// The index of the face in the compressed grid
const int faceIdx = *cellFaceIter;
// the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
switch(faceTag)
{
case 0: faceDirection = Opm::FaceDir::XMinus;
break;
case 1: faceDirection = Opm::FaceDir::XPlus;
break;
case 2: faceDirection = Opm::FaceDir::YMinus;
break;
case 3: faceDirection = Opm::FaceDir::YPlus;
break;
case 4: faceDirection = Opm::FaceDir::ZMinus;
break;
case 5: faceDirection = Opm::FaceDir::ZPlus;
break;
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Fetkovich problem. Make sure faceTag is correctly defined");
}
if (faceDirection == connection.reservoir_face_dir.at(idx))
{
faceArea_connected_.at(idx) = getFaceArea(faceCells, ugrid, faceIdx, idx, connection);
denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) );
}
}
auto cellCenter = grid.getCellCenter(cell_idx_.at(idx));
cell_depth_.at(idx) = cellCenter[2];
}
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
{
alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
0.
: ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas;
}
}
inline void calculateAquiferCondition()
{
int pvttableIdx = aqufetp_data_.pvttableID - 1;
rhow_.resize(cell_idx_.size(),0.);
if (!aqufetp_data_.p0)
{
pa0_ = calculateReservoirEquilibrium();
}
else
{
pa0_ = *(aqufetp_data_.p0);
}
aquifer_pressure_ = pa0_ ;
// use the thermodynamic state of the first active cell as a
// reference. there might be better ways to do this...
ElementContext elemCtx(ebos_simulator_);
auto elemIt = 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 = pa0_;
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
mu_w_ = mu_w_aquifer.value();
}
inline Scalar calculateReservoirEquilibrium()
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(ebos_simulator_);
const auto& gridView = ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back( (water_pressure_reservoir - rhow_[idx].value()*gravity_()*(cell_depth_[idx] - aqufetp_data_.d0))*alphai_[idx] );
}
// 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();
return aquifer_pres_avg;
}
}; //Class AquiferFetkovich
} // namespace Opm
#endif

226
opm/autodiff/AquiferInterface.hpp Normal file
View File

@ -0,0 +1,226 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifetp.hpp>
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp>
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <vector>
#include <algorithm>
#include <unordered_map>
namespace Opm
{
template<typename TypeTag>
class AquiferInterface
{
public:
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
static const int numEq = BlackoilIndices::numEq;
typedef double Scalar;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef Opm::BlackOilFluidState<Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, BlackoilIndices::numPhases> FluidState;
static const auto waterCompIdx = FluidSystem::waterCompIdx;
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
// Constructor
AquiferInterface( const Aquancon::AquanconOutput& connection,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator)
: connection_(connection)
, ebos_simulator_(ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed)
{}
// Deconstructor
virtual ~AquiferInterface() {}
void initialSolutionApplied()
{
initQuantities(connection_);
}
void beginTimeStep()
{
ElementContext elemCtx(ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin<0>();
const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
int cellIdx = elemCtx.globalSpaceIndex(0, 0);
int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updateIntensiveQuantities(0);
const auto& iq = elemCtx.intensiveQuantities(0, 0);
pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx));
}
}
template <class Context>
void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx)
{
unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
return;
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to
// IntensiveQuantities of that particular cell_id
const IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
// This is the pressure at td + dt
updateCellPressure(pressure_current_,idx,intQuants);
updateCellDensity(idx,intQuants);
calculateInflowRate(idx, context.simulator());
rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] +=
Qai_[idx]/context.dofVolume(spaceIdx, timeIdx);
}
protected:
inline Scalar gravity_() const
{
return ebos_simulator_.problem().gravity()[2];
}
inline void initQuantities(const Aquancon::AquanconOutput& connection)
{
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0
W_flux_ = 0.;
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections function
initializeConnections(connection);
calculateAquiferCondition();
calculateAquiferConstants();
pressure_previous_.resize(cell_idx_.size(), 0.);
pressure_current_.resize(cell_idx_.size(), 0.);
Qai_.resize(cell_idx_.size(), 0.0);
}
inline void updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
}
inline void updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value();
}
inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
rhow_.at(idx) = fs.density(waterPhaseIdx);
}
template<class faceCellType, class ugridType>
inline const double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid,
const int faceIdx, const int idx,
const Aquancon::AquanconOutput& connection) const
{
// 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
// face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining)
double faceArea = 0.;
const auto cellNeighbour0 = faceCells(faceIdx,0);
const auto cellNeighbour1 = faceCells(faceIdx,1);
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))?
defaultFaceArea :
*(connection.influx_coeff.at(idx));
faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea;
if (cellNeighbour1 == 0){
faceArea = (cellNeighbour0 < 0)? faceArea : 0.;
}
else if (cellNeighbour0 == 0){
faceArea = (cellNeighbour1 < 0)? faceArea : 0.;
}
return faceArea;
}
virtual void endTimeStep() = 0;
const Aquancon::AquanconOutput connection_;
const Simulator& ebos_simulator_;
const std::unordered_map<int, int> cartesian_to_compressed_;
// Grid variables
std::vector<size_t> cell_idx_;
std::vector<Scalar> faceArea_connected_;
std::vector<int> cellToConnectionIdx_;
// Quantities at each grid id
std::vector<Scalar> cell_depth_;
std::vector<Scalar> pressure_previous_;
std::vector<Eval> pressure_current_;
std::vector<Eval> Qai_;
std::vector<Eval> rhow_;
std::vector<Scalar> alphai_;
Scalar mu_w_; //water viscosity
Scalar Tc_; // Time constant
Scalar pa0_; // initial aquifer pressure
Eval W_flux_;
virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0;
virtual Scalar dpai(int idx) = 0;
virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
virtual void calculateAquiferCondition() = 0;
virtual void calculateAquiferConstants() = 0;
virtual Scalar calculateReservoirEquilibrium() =0;
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
};
} // namespace Opm
#endif

4
opm/autodiff/BlackoilAquiferModel_impl.hpp
View File

@ -134,7 +134,7 @@ namespace Opm {
for (size_t i = 0; i < aquifersData.size(); ++i) for (size_t i = 0; i < aquifersData.size(); ++i)
{ {
aquifers_CarterTracy.push_back( aquifers_CarterTracy.push_back(
AquiferCarterTracy<TypeTag> (aquifersData.at(i), aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_) AquiferCarterTracy<TypeTag> (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i))
); );
} }
} }
@ -161,7 +161,7 @@ namespace Opm {
for (size_t i = 0; i < aquifersData.size(); ++i) for (size_t i = 0; i < aquifersData.size(); ++i)
{ {
aquifers_Fetkovich.push_back( aquifers_Fetkovich.push_back(
AquiferFetkovich<TypeTag> (aquifersData.at(i), aquifer_connection.at(i),cartesian_to_compressed_, this->simulator_) AquiferFetkovich<TypeTag> (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i))
); );
} }
} }