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

Runs first iterations of two-phase case plausibly

Browse Source
This commit is contained in:
babrodtk 2016-04-21 16:19:45 +02:00
parent 2ac514a1ea
commit a1504a2bff
5 changed files with 103 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/AutoDiffBlock.hpp
View File

@ -292,6 +292,7 @@ namespace Opm
/// Elementwise operator -
AutoDiffBlock operator-(const AutoDiffBlock& rhs) const
{
/*
assert(val_.rows() > 0 || rhs.val_.rows() > 0);
if (val_.rows() == 0) {
return rhs*(-1.0);
@ -299,6 +300,7 @@ namespace Opm
if (rhs.val_.rows() == 0) {
return *this;
}
*/
if (jac_.empty() && rhs.jac_.empty()) {
return constant(val_ - rhs.val_);

107
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -56,6 +56,7 @@
#include <iomanip>
#include <limits>
#include <vector>
#include <algorithm>
//#include <fstream>
// A debugging utility.
@ -184,10 +185,21 @@ namespace detail {
{ 1.1169, 1.0031, 0.0031 }, // the default magic numbers
false } )
, terminal_output_ (terminal_output)
, material_name_{ "Water", "Oil", "Gas" }
, material_name_(0)
, current_relaxation_(1.0)
{
assert(numMaterials() == 3); // Due to the material_name_ init above.
if (active_[Water]) {
material_name_.push_back("Water");
}
if (active_[Oil]) {
material_name_.push_back("Oil");
}
if (active_[Gas]) {
material_name_.push_back("Gas");
}
assert(numMaterials() == std::accumulate(active_.begin(), active_.end(), 0)); // Due to the material_name_ init above.
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
@ -662,7 +674,7 @@ namespace detail {
{
const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ]
: ADB::constant(V::Zero(nc, 1), sg.blockPattern()));
: ADB::null());
state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
}
@ -686,8 +698,8 @@ namespace detail {
// Compute phase pressures also if gas phase is not active
const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ]
: ADB::constant(V::Zero(nc, 1), so.blockPattern()));
const ADB& sg = ADB::constant(V::Zero(nc, 1), so.blockPattern());
: ADB::null());
const ADB& sg = ADB::null();
state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
}
@ -1590,37 +1602,43 @@ namespace detail {
sg = sg_old - step * dsg;
}
assert(active_[Oil]);
const int pos = pu.phase_pos[ Oil ];
const V so_old = s_old.col(pos);
so = so_old - step * dso;
}
// Appleyard chop process.
auto ixg = sg < 0;
for (int c = 0; c < nc; ++c) {
if (ixg[c]) {
sw[c] = sw[c] / (1-sg[c]);
so[c] = so[c] / (1-sg[c]);
sg[c] = 0;
if (active_[Gas]) {
auto ixg = sg < 0;
for (int c = 0; c < nc; ++c) {
if (ixg[c]) {
sw[c] = sw[c] / (1-sg[c]);
so[c] = so[c] / (1-sg[c]);
sg[c] = 0;
}
}
}
auto ixo = so < 0;
for (int c = 0; c < nc; ++c) {
if (ixo[c]) {
sw[c] = sw[c] / (1-so[c]);
sg[c] = sg[c] / (1-so[c]);
so[c] = 0;
if (active_[Oil]) {
auto ixo = so < 0;
for (int c = 0; c < nc; ++c) {
if (ixo[c]) {
sw[c] = sw[c] / (1-so[c]);
sg[c] = sg[c] / (1-so[c]);
so[c] = 0;
}
}
}
auto ixw = sw < 0;
for (int c = 0; c < nc; ++c) {
if (ixw[c]) {
so[c] = so[c] / (1-sw[c]);
sg[c] = sg[c] / (1-sw[c]);
sw[c] = 0;
if (active_[Water]) {
auto ixw = sw < 0;
for (int c = 0; c < nc; ++c) {
if (ixw[c]) {
so[c] = so[c] / (1-sw[c]);
sg[c] = sg[c] / (1-sw[c]);
sw[c] = 0;
}
}
}
@ -1630,18 +1648,21 @@ namespace detail {
//sg = sg / sumSat;
// Update the reservoir_state
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
if (active_[Water]) {
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
}
}
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
if (active_[Gas]) {
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
}
}
if (active_[ Oil ]) {
const int pos = pu.phase_pos[ Oil ];
for (int c = 0; c < nc; ++c) {
reservoir_state.saturation()[c*np + pos] = so[c];
reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
}
}
@ -1887,7 +1908,9 @@ namespace detail {
// The reference pressure is always the liquid phase (oil) pressure.
if (phaseIdx == BlackoilPhases::Liquid)
continue;
pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
if (active_[phaseIdx]) {
pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
}
}
// Since pcow = po - pw, but pcog = pg - po,
@ -1896,10 +1919,12 @@ namespace detail {
// pg = po + pcgo
// This is an unfortunate inconsistency, but a convention we must handle.
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
if (phaseIdx == BlackoilPhases::Aqua) {
pressure[phaseIdx] = po - pressure[phaseIdx];
} else {
pressure[phaseIdx] += po;
if (active_[phaseIdx]) {
if (phaseIdx == BlackoilPhases::Aqua) {
pressure[phaseIdx] = po - pressure[phaseIdx];
} else {
pressure[phaseIdx] += po;
}
}
}
@ -2660,13 +2685,17 @@ namespace detail {
BlackoilModelBase<Grid, WellModel, Implementation>::
updatePhaseCondFromPrimalVariable()
{
if (! active_[Gas]) {
OPM_THROW(std::logic_error, "updatePhaseCondFromPrimarVariable() logic requires active gas phase.");
}
const int nc = primalVariable_.size();
const int nc = Opm::AutoDiffGrid::numCells(grid_);
isRs_ = V::Zero(nc);
isRv_ = V::Zero(nc);
isSg_ = V::Zero(nc);
if (! (active_[Gas] && active_[Oil])) {
// updatePhaseCondFromPrimarVariable() logic requires active gas and oil phase.
phaseCondition_.assign(nc, PhasePresence());
return;
}
for (int c = 0; c < nc; ++c) {
phaseCondition_[c] = PhasePresence(); // No free phases.
phaseCondition_[c].setFreeWater(); // Not necessary for property calculation usage.

20
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -550,18 +550,19 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
pLad.value = po.value()[i];
TLad.value = T.value()[i];
if (phase_usage_.phase_used[Gas]) {
//RS/RV only makes sense when gas phase is active
if (cond[i].hasFreeGas()) {
bLad = oilPvt_->saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
//RS/RV only makes sense when gas phase is active
if (cond[i].hasFreeGas()) {
bLad = oilPvt_->saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
}
else {
if (rs.size() == 0) {
// If FIXME
RsLad.value = 0.0;
}
else {
RsLad.value = rs.value()[i];
bLad = oilPvt_->inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
}
}
else {
bLad = 0.0;
bLad = oilPvt_->inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
}
b[i] = bLad.value;
@ -867,9 +868,6 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
std::vector<ADB> adbCapPressures;
adbCapPressures.reserve(3);
const ADB* s[3] = { &sw, &so, &sg };
//const std::vector<int>& block_pattern = sw.blockPattern();
assert(sw.blockPattern() == so.blockPattern());
assert(sw.blockPattern() == sg.blockPattern());
for (int phase1 = 0; phase1 < 3; ++phase1) {
if (phase_usage_.phase_used[phase1]) {
const int phase1_pos = phase_usage_.phase_pos[phase1];

16
opm/autodiff/StandardWells_impl.hpp
View File

@ -219,18 +219,20 @@ namespace Opm
}
assert(active[Oil]);
const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rs = subset(state.rs, well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null();
const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
// const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null();
const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
}
if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
const ADB perf_rv = subset(state.rv, well_cells);
const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
// const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
}

25
opm/simulators/thresholdPressures.hpp
View File

@ -110,20 +110,25 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
// the inverse formation volume factors
std::vector<PhasePresence> cond(numCells);
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
if (pu.phase_used[BlackoilPhases::Aqua] && sw > 0.0) {
cond[cellIdx].setFreeWater();
if (pu.phase_used[BlackoilPhases::Aqua]) {
const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
if (sw > 0.0) {
cond[cellIdx].setFreeWater();
}
}
if (pu.phase_used[BlackoilPhases::Liquid] && so > 0.0) {
cond[cellIdx].setFreeOil();
if (pu.phase_used[BlackoilPhases::Liquid]) {
const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
if (so > 0.0) {
cond[cellIdx].setFreeOil();
}
}
if (pu.phase_used[BlackoilPhases::Vapour] && sg > 0.0) {
cond[cellIdx].setFreeGas();
if (pu.phase_used[BlackoilPhases::Vapour]) {
const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
if (sg > 0.0) {
cond[cellIdx].setFreeGas();
}
}
}