Handle combination of threshold_pressure and nnc

The nncs threshold pressures are calculated and appended to the face
threshold pressures
This commit is contained in:
Tor Harald Sandve
2015-12-08 10:56:42 +01:00
parent 99ddc46318
commit e3393c5ee9
8 changed files with 26 additions and 11 deletions

View File

@@ -413,7 +413,8 @@ try
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
SimulatorFullyImplicitBlackoil< Grid > simulator(param,
grid,
geoprops,

View File

@@ -413,6 +413,8 @@ try
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
SimulatorFullyImplicitBlackoilMultiSegment< Grid > simulator(param,
grid,

View File

@@ -287,6 +287,8 @@ try
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, *grid->c_grid(), state, *props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid(), maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geology.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
Opm::BlackoilOutputWriter
outputWriter(cGrid, param, eclipseState, pu,

View File

@@ -383,6 +383,8 @@ try
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
SimulatorFullyImplicitBlackoilSolvent< Grid > simulator(param,
grid,

View File

@@ -298,7 +298,7 @@ namespace Opm {
ModelParameters param_;
bool use_threshold_pressure_;
bool wells_active_;
V threshold_pressures_by_interior_face_;
V threshold_pressures_by_connection_;
std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_;

View File

@@ -370,19 +370,27 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::
setThresholdPressures(const std::vector<double>& threshold_pressures_by_face)
setThresholdPressures(const std::vector<double>& threshold_pressures)
{
const int num_faces = AutoDiffGrid::numFaces(grid_);
if (int(threshold_pressures_by_face.size()) != num_faces) {
OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces.");
const int num_nnc = geo_.nnc().numNNC();
const int num_connections = num_faces + num_nnc;
if (int(threshold_pressures.size()) != num_connections) {
OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures input, must be equal to number of faces + nncs.");
}
use_threshold_pressure_ = true;
// Map to interior faces.
const int num_ifaces = ops_.internal_faces.size();
threshold_pressures_by_interior_face_.resize(num_ifaces);
threshold_pressures_by_connection_.resize(num_ifaces + num_nnc);
for (int ii = 0; ii < num_ifaces; ++ii) {
threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]];
threshold_pressures_by_connection_[ii] = threshold_pressures[ops_.internal_faces[ii]];
}
// Handle NNCs
// Note: the nnc threshold pressures is appended after the face threshold pressures
for (int ii = 0; ii < num_nnc; ++ii) {
threshold_pressures_by_connection_[ii + num_ifaces] = threshold_pressures[ii + num_faces];
}
}
@@ -2421,14 +2429,14 @@ namespace detail {
// threshold, that shall have zero flow. Storing the bool
// Array as a V (a double Array) with 1 and 0 elements, a
// 1 where flow is allowed, a 0 where it is not.
const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast<double>();
const V high_potential = (dp.value().abs() >= threshold_pressures_by_connection_).template cast<double>();
// Create a sparse vector that nullifies the low potential elements.
const M keep_high_potential(high_potential.matrix().asDiagonal());
// Find the current sign for the threshold modification
const V sign_dp = sign(dp.value());
const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_;
const V threshold_modification = sign_dp * threshold_pressures_by_connection_;
// Modify potential and nullify where appropriate.
dp = keep_high_potential * (dp - threshold_modification);

View File

@@ -123,7 +123,7 @@ namespace Opm {
using Base::has_vapoil_;
using Base::param_;
using Base::use_threshold_pressure_;
using Base::threshold_pressures_by_interior_face_;
using Base::threshold_pressures_by_connection_;
using Base::rq_;
using Base::phaseCondition_;
using Base::well_perforation_pressure_diffs_;

View File

@@ -169,7 +169,7 @@ namespace Opm {
using Base::has_vapoil_;
using Base::param_;
using Base::use_threshold_pressure_;
using Base::threshold_pressures_by_interior_face_;
using Base::threshold_pressures_by_connection_;
using Base::rq_;
using Base::phaseCondition_;
using Base::well_perforation_pressure_diffs_;