Add applyThresholdPressures() method and usage.

Note that this commit does not introduce any way to set
use_threshold_pressure to true, so the new code is not run.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-08-27 10:27:49 +02:00
parent 102881bff8
commit f7fa3488cb
2 changed files with 42 additions and 1 deletions

View File

@ -155,6 +155,8 @@ namespace Opm {
double relax_increment_;
double relax_rel_tol_;
int max_iter_;
bool use_threshold_pressure_;
V threshold_pressures_by_face_;
std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_;
@ -223,6 +225,8 @@ namespace Opm {
const ADB& p ,
const SolutionState& state );
void applyThresholdPressures(ADB& dp);
double
residualNorm() const;

View File

@ -190,6 +190,7 @@ namespace {
, relax_increment_ (0.1)
, relax_rel_tol_ (0.2)
, max_iter_ (15)
, use_threshold_pressure_(false)
, rq_ (fluid.numPhases())
, phaseCondition_(AutoDiffGrid::numCells(grid))
, residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
@ -1577,7 +1578,11 @@ namespace {
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
const ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(dp);
}
head = transi*dp;
//head = transi*(ops_.ngrad * phasePressure) + gflux;
@ -1595,6 +1600,38 @@ namespace {
template<class T>
void
FullyImplicitBlackoilSolver<T>::applyThresholdPressures(ADB& dp)
{
// We support reversible threshold pressures only.
// Method: if the potential difference is lower (in absolute
// value) than the threshold for any face, then the potential
// (and derivatives) is set to zero. If it is above the
// threshold, the threshold pressure is subtracted from the
// absolute potential (the potential is moved towards zero).
// Identify the set of faces where the potential is under the
// threshold, that shall have zero flow. Storing the bool
// Array as a V (a double Array) with 1 and 0 elements.
const V low_potential = (dp.value().abs() < threshold_pressures_by_face_).template cast<double>();
// Create a sparse vector that nullifies the low potential elements.
const M nullify_low_potential = spdiag(low_potential);
// The threshold modification must have the sign that reduces the
// absolute value of the potential.
const V sign_for_mod = (dp.value() < 0).template cast<double>();
const V threshold_modification = sign_for_mod * threshold_pressures_by_face_;
// Modify potential and nullify where appropriate.
dp = nullify_low_potential * (dp + threshold_modification);
}
template<class T>
double
FullyImplicitBlackoilSolver<T>::residualNorm() const