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

Base Class and AquiferCarterTracy and AquiferFetkovich implementation

Browse Source 
Reordering of Aquifer Codes
This commit is contained in:
WesselDeZeeuw 2019-01-07 11:44:33 +01:00
parent 863dcf1328
commit a72d61cb50
4 changed files with 519 additions and 602 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:
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AquiferInterface.hpp
opm/autodiff/AquiferCarterTracy.hpp
opm/autodiff/AquiferFetkovich.hpp
opm/autodiff/BlackoilAmg.hpp
opm/autodiff/BlackoilDetails.hpp
opm/autodiff/BlackoilModelParametersEbos.hpp

270
opm/autodiff/AquiferCarterTracy.hpp
View File

@ -21,141 +21,48 @@
#ifndef OPM_AQUIFERCT_HEADER_INCLUDED
#define OPM_AQUIFERCT_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/AquiferCT.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>
#include <opm/autodiff/AquiferInterface.hpp>
namespace Opm
{
template<typename TypeTag>
class AquiferCarterTracy
class AquiferCarterTracy: public AquiferInterface<TypeTag>
{
public:
typedef AquiferInterface<TypeTag> Base;
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;
using typename Base::Simulator;
using typename Base::ElementContext;
using typename Base::FluidSystem;
using typename Base::BlackoilIndices;
using typename Base::RateVector;
using typename Base::IntensiveQuantities;
using typename Base::Eval;
using typename Base::Scalar;
using typename Base::FluidState;
using Base::waterCompIdx;
using Base::waterPhaseIdx;
AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data,
const Aquancon::AquanconOutput& connection,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator)
: ebos_simulator_ (ebosSimulator)
: Base(connection, cartesian_to_compressed, ebosSimulator)
, aquct_data_(aquct_data)
, cartesian_to_compressed_(cartesian_to_compressed)
, 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()
{
for (const auto& Qai: Qai_) {
W_flux_ += Qai*ebos_simulator_.timeStepSize();
for (const auto& Qai: Base::Qai_) {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize();
}
}
private:
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_;
protected:
// Variables constants
const AquiferCT::AQUCT_data aquct_data_;
Scalar mu_w_; //water viscosity
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)
{
@ -164,57 +71,22 @@ namespace Opm
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);
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()) / Tc_;
const Scalar td = simulator.time() / Tc_;
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/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime );
b = beta_ / (Tc_ * ( PItd - td*PItdprime));
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
@ -222,7 +94,7 @@ namespace Opm
{
Scalar a, b;
calculateEqnConstants(a,b,idx,simulator);
Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) );
Base::Qai_.at(idx) = Base::alphai_.at(idx)*( a - b * ( Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx) ) );
}
inline void calculateAquiferConstants()
@ -233,56 +105,31 @@ namespace Opm
* aquct_data_.C_t
* aquct_data_.r_o * aquct_data_.r_o;
// We calculate the time constant
Tc_ = mu_w_ * aquct_data_.phi_aq
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 );
}
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
inline void initializeConnections(const Aquancon::AquanconOutput& connection)
{
const auto& eclState = ebos_simulator_.vanguard().eclState();
const auto& ugrid = ebos_simulator_.vanguard().grid();
const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
cell_idx_ = connection.global_index;
Base::cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell();
assert( cell_idx_ == connection.global_index);
assert( (cell_idx_.size() <= connection.influx_coeff.size()) );
assert( Base::cell_idx_ == connection.global_index);
assert( (Base::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(), aquct_data_.d0);
alphai_.resize(cell_idx_.size(), 1.0);
faceArea_connected_.resize(cell_idx_.size(),0.0);
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);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
@ -292,11 +139,11 @@ namespace Opm
// 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)
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 = cartesian_to_compressed_.at(cell_idx_[idx]);
cellToConnectionIdx_[cell_index] = 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];
for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
@ -326,21 +173,21 @@ namespace Opm
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) );
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(cell_idx_.at(idx));
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());
for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
for (size_t idx = 0; idx < Base::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.
: ( 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;
}
}
@ -349,21 +196,21 @@ namespace Opm
int pvttableIdx = aquct_data_.pvttableID - 1;
rhow_.resize(cell_idx_.size(),0.);
Base::rhow_.resize(Base::cell_idx_.size(),0.);
if (!aquct_data_.p0)
{
pa0_ = calculateReservoirEquilibrium();
Base::pa0_ = calculateReservoirEquilibrium();
}
else
{
pa0_ = *(aquct_data_.p0);
Base::pa0_ = *(aquct_data_.p0);
}
// 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>();
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);
@ -375,11 +222,11 @@ namespace Opm
fs_aquifer.assign( iq0.fluidState() );
Eval temperature_aq, pa0_mean;
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);
mu_w_ = mu_w_aquifer.value();
Base::mu_w_ = mu_w_aquifer.value();
}
@ -390,8 +237,8 @@ namespace Opm
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(ebos_simulator_);
const auto& gridView = ebos_simulator_.gridView();
ElementContext elemCtx(Base::ebos_simulator_);
const auto& gridView = Base::ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@ -399,7 +246,7 @@ namespace Opm
elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = cellToConnectionIdx_[cellIdx];
int idx = Base::cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
@ -408,20 +255,15 @@ namespace Opm
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] - aquct_data_.d0))*alphai_[idx] );
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] );
}
// 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;
}
const Aquancon::AquanconOutput connection_;
std::vector<int> cellToConnectionIdx_;
}; // class AquiferCarterTracy
} // namespace Opm
#endif

305
opm/autodiff/AquiferFetkovich.hpp
View File

@ -21,240 +21,70 @@ along with OPM. If not, see <http://www.gnu.org/licenses/>.
#ifndef OPM_AQUIFETP_HEADER_INCLUDED
#define OPM_AQUIFETP_HEADER_INCLUDED
#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>
#include <opm/autodiff/AquiferInterface.hpp>
namespace Opm
{
template<typename TypeTag>
class AquiferFetkovich
class AquiferFetkovich: public AquiferInterface<TypeTag>
{
public:
typedef AquiferInterface<TypeTag> Base;
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) };
using typename Base::Simulator;
using typename Base::ElementContext;
using typename Base::FluidSystem;
using typename Base::BlackoilIndices;
using typename Base::RateVector;
using typename Base::IntensiveQuantities;
using typename Base::Eval;
using typename Base::Scalar;
using typename Base::FluidState;
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;
using Base::waterCompIdx;
using Base::waterPhaseIdx;
AquiferFetkovich( const Aquifetp::AQUFETP_data& aqufetp_data,
const Aquancon::AquanconOutput& connection,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator)
: ebos_simulator_ (ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed)
: Base(connection, cartesian_to_compressed, ebosSimulator)
, aqufetp_data_(aqufetp_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()
{
for (const auto& Qai: Qai_) {
W_flux_ += Qai*ebos_simulator_.timeStepSize();
for (const auto& Qai: Base::Qai_) {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize();
aquifer_pressure_ = aquiferPressure();
}
}
private:
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_;
// 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_;
std::vector<int> cellToConnectionIdx_;
// Variables constants
protected:
// Aquifer Fetkovich Specific Variables
const Aquifetp::AQUFETP_data aqufetp_data_;
const Aquancon::AquanconOutput connection_;
Scalar aquifer_pressure_; // aquifer
Scalar mu_w_; //water viscosity
Scalar Tc_; // Time Constant
Scalar pa0_; // initial aquifer pressure
Scalar aquifer_pressure_; // aquifer pressure
Eval W_flux_;
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();
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 = aquifer_pressure_ + rhow_.at(idx).value()*gravity_()*(cell_depth_.at(idx) - aqufetp_data_.d0) - pressure_current_.at(idx).value() ;
return dp;
}
// This function implements Eq 5.12 of the EclipseTechnicalDescription
inline Scalar aquiferPressure()
{
Scalar Flux = W_flux_.value();
Scalar pa_ = pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 );
return pa_;
}
// This function implements Eq 5.14 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator)
{
Tc_ = ( aqufetp_data_.C_t * aqufetp_data_.V0 ) / aqufetp_data_.J ;
Scalar td_Tc_ = simulator.timeStepSize() / Tc_ ;
Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_;
Qai_.at(idx) = alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_;
}
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
inline void initializeConnections(const Aquancon::AquanconOutput& connection)
{
const auto& eclState = ebos_simulator_.vanguard().eclState();
const auto& ugrid = ebos_simulator_.vanguard().grid();
const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
cell_idx_ = connection.global_index;
Base::cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell();
assert( cell_idx_ == connection.global_index);
assert( (cell_idx_.size() == connection.influx_coeff.size()) );
assert( Base::cell_idx_ == connection.global_index);
assert( (Base::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);
Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0);
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);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
@ -264,11 +94,11 @@ namespace Opm
// 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)
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 = cartesian_to_compressed_.at(cell_idx_[idx]);
cellToConnectionIdx_[cell_index] = 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];
for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
{
@ -292,49 +122,72 @@ namespace Opm
break;
case 5: faceDirection = Opm::FaceDir::ZPlus;
break;
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Fetkovich problem. Make sure faceTag is correctly defined");
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer 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) );
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(cell_idx_.at(idx));
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());
for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
for (size_t idx = 0; idx < Base::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.
: ( 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 Scalar dpai(int idx)
{
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() ;
return dp;
}
// This function implements Eq 5.12 of the EclipseTechnicalDescription
inline Scalar aquiferPressure()
{
Scalar Flux = Base::W_flux_.value();
Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 );
return pa_;
}
// This function implements Eq 5.14 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator)
{
Base::Tc_ = ( aqufetp_data_.C_t * aqufetp_data_.V0 ) / aqufetp_data_.J ;
Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ;
Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_;
Base::Qai_.at(idx) = Base::alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_;
}
inline void calculateAquiferCondition()
{
int pvttableIdx = aqufetp_data_.pvttableID - 1;
rhow_.resize(cell_idx_.size(),0.);
Base::rhow_.resize(Base::cell_idx_.size(),0.);
if (!aqufetp_data_.p0)
{
pa0_ = calculateReservoirEquilibrium();
Base::pa0_ = calculateReservoirEquilibrium();
}
else
{
pa0_ = *(aqufetp_data_.p0);
Base::pa0_ = *(aqufetp_data_.p0);
}
aquifer_pressure_ = pa0_ ;
aquifer_pressure_ = Base::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>();
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
@ -342,11 +195,9 @@ namespace Opm
fs_aquifer.assign( iq0.fluidState() );
Eval temperature_aq, pa0_mean;
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);
mu_w_ = mu_w_aquifer.value();
Base::mu_w_ = mu_w_aquifer.value();
}
inline Scalar calculateReservoirEquilibrium()
@ -355,16 +206,15 @@ namespace Opm
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(ebos_simulator_);
const auto& gridView = ebos_simulator_.gridView();
ElementContext elemCtx(Base::ebos_simulator_);
const auto& gridView = Base::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];
int idx = Base::cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
@ -373,8 +223,8 @@ namespace Opm
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] );
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] );
}
// We take the average of the calculated equilibrium pressures.
@ -383,5 +233,4 @@ namespace Opm
}
}; //Class AquiferFetkovich
} // namespace Opm
#endif

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

@ -0,0 +1,224 @@
/*
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)
: ebos_simulator_(ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed)
, connection_(connection)
{}
// 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);
}
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();
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;
protected:
const Simulator& ebos_simulator_;
const Aquancon::AquanconOutput connection_;
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 Scalar calculateReservoirEquilibrium() =0;
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
};
} // namespace Opm
#endif