mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
add capPress functionality for PEDs, just use phase pressure to compute
laplace term, all the properties are computed by reference pressure which maybe oil pressure in OPM.
This commit is contained in:
Liu Ming
parent
0055ae1def
commit
c4d567c5e4
@ -584,5 +584,64 @@ namespace Opm
|
||||
return relperms;
|
||||
}
|
||||
|
||||
std::vector<ADB> BlackoilPropsAd::capPress(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const
|
||||
|
||||
{
|
||||
const int numCells = cells.size();
|
||||
const int numActivePhases = numPhases();
|
||||
const int numBlocks = so.numBlocks();
|
||||
|
||||
Block activeSat(numCells, numActivePhases);
|
||||
if (pu_.phase_used[Water]) {
|
||||
assert(sw.value().size() == numCells);
|
||||
activeSat.col(pu_.phase_pos[Water]) = sw.value();
|
||||
}
|
||||
if (pu_.phase_used[Oil]) {
|
||||
assert(so.value().size() == numCells);
|
||||
activeSat.col(pu_.phase_pos[Oil]) = so.value();
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
||||
}
|
||||
if (pu_.phase_used[Gas]) {
|
||||
assert(sg.value().size() == numCells);
|
||||
activeSat.col(pu_.phase_pos[Gas]) = sg.value();
|
||||
}
|
||||
|
||||
Block pc(numCells, numActivePhases);
|
||||
Block dpc(numCells, numActivePhases*numActivePhases);
|
||||
props_.capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
|
||||
|
||||
std::vector<ADB> adbCapPressures;
|
||||
adbCapPressures.reserve(3);
|
||||
const ADB* s[3] = { &sw, &so, &sg };
|
||||
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
||||
if (pu_.phase_used[phase1]) {
|
||||
const int phase1_pos = pu_.phase_pos[phase1];
|
||||
std::vector<ADB::M> jacs(numBlocks);
|
||||
for (int block = 0; block < numBlocks; ++block) {
|
||||
jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
|
||||
}
|
||||
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
||||
if (!pu_.phase_used[phase2])
|
||||
continue;
|
||||
const int phase2_pos = pu_.phase_pos[phase2];
|
||||
// Assemble dpc1/ds2.
|
||||
const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
||||
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
|
||||
for (int block = 0; block < numBlocks; ++block) {
|
||||
jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
|
||||
}
|
||||
}
|
||||
adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
|
||||
} else {
|
||||
adbCapPressures.emplace_back(ADB::null());
|
||||
}
|
||||
}
|
||||
return adbCapPressures;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
@ -237,6 +237,18 @@ namespace Opm
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
/// Capillary pressure for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n capillary pressure values,
|
||||
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
|
||||
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
|
||||
std::vector<ADB> capPress(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
|
||||
@ -678,5 +678,63 @@ namespace Opm
|
||||
return relperms;
|
||||
}
|
||||
|
||||
std::vector<ADB> BlackoilPropsAdFromDeck::capPress(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int numCells = cells.size();
|
||||
const int numActivePhases = numPhases();
|
||||
const int numBlocks = so.numBlocks();
|
||||
|
||||
Block activeSat(numCells, numActivePhases);
|
||||
if (phase_usage_.phase_used[Water]) {
|
||||
assert(sw.value().size() == numCells);
|
||||
activeSat.col(phase_usage_.phase_pos[Water]) = sw.value();
|
||||
}
|
||||
if (phase_usage_.phase_used[Oil]) {
|
||||
assert(so.value().size() == numCells);
|
||||
activeSat.col(phase_usage_.phase_pos[Oil]) = so.value();
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
||||
}
|
||||
if (phase_usage_.phase_used[Gas]) {
|
||||
assert(sg.value().size() == numCells);
|
||||
activeSat.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
||||
}
|
||||
|
||||
Block pc(numCells, numActivePhases);
|
||||
Block dpc(numCells, numActivePhases*numActivePhases);
|
||||
satprops_->capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
|
||||
|
||||
std::vector<ADB> adbCapPressures;
|
||||
adbCapPressures.reserve(3);
|
||||
const ADB* s[3] = { &sw, &so, &sg };
|
||||
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
||||
if (phase_usage_.phase_used[phase1]) {
|
||||
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
||||
std::vector<ADB::M> jacs(numBlocks);
|
||||
for (int block = 0; block < numBlocks; ++block) {
|
||||
jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
|
||||
}
|
||||
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
||||
if (!phase_usage_.phase_used[phase2])
|
||||
continue;
|
||||
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
||||
// Assemble dpc1/ds2.
|
||||
const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
||||
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
|
||||
for (int block = 0; block < numBlocks; ++block) {
|
||||
jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
|
||||
}
|
||||
}
|
||||
adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
|
||||
} else {
|
||||
adbCapPressures.emplace_back(ADB::null());
|
||||
}
|
||||
}
|
||||
return adbCapPressures;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
@ -238,6 +238,19 @@ namespace Opm
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
/// Capillary pressure for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n capillary pressure values,
|
||||
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
|
||||
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
|
||||
std::vector<ADB> capPress(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
private:
|
||||
RockFromDeck rock_;
|
||||
|
||||
@ -243,6 +243,21 @@ namespace Opm
|
||||
const ADB& sg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
|
||||
/// Capillary pressure for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n capillary pressure values,
|
||||
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
|
||||
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
|
||||
virtual
|
||||
std::vector<ADB> capPress(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
@ -504,7 +504,6 @@ namespace {
|
||||
// for each active phase.
|
||||
const V trans = subset(geo_.transmissibility(), ops_.internal_faces);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
// const ADB cmax = computeCmax(state.concentration);
|
||||
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
|
||||
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
|
||||
const ADB mc = computeMc(state);
|
||||
@ -815,6 +814,33 @@ namespace {
|
||||
|
||||
|
||||
|
||||
std::vector<ADB>
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
computePressures(const SolutionState& state) const
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const std::vector<int>& bpat = state.pressure.blockPattern();
|
||||
|
||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||
|
||||
const ADB sw = state.saturation[0];
|
||||
|
||||
const ADB so = state.saturation[1];
|
||||
|
||||
const ADB sg = null;
|
||||
|
||||
// convert the pressure offsets to the capillary pressures
|
||||
std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
|
||||
pressure[0] = pressure[0] - pressure[0];
|
||||
|
||||
// add the total pressure to the capillary pressures
|
||||
for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) {
|
||||
pressure[phaseIdx] += state.pressure;
|
||||
}
|
||||
|
||||
return pressure;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void
|
||||
@ -833,12 +859,13 @@ namespace {
|
||||
rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
|
||||
const ADB mu_o = fluidViscosity(1, state.pressure, cells_);
|
||||
rq_[1].mob = tr_mult * kro / mu_o;
|
||||
std::vector<ADB> press = computePressures(state);
|
||||
for (int phase = 0; phase < 2; ++phase) {
|
||||
const ADB rho = fluidDensity(phase, state.pressure, cells_);
|
||||
ADB& head = rq_[ phase ].head;
|
||||
// compute gravity potensial using the face average as in eclipse and MRST
|
||||
const ADB rhoavg = ops_.caver * rho;
|
||||
const ADB dp = ops_.ngrad * state.pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
const ADB dp = ops_.ngrad * press[phase] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
head = transi*dp;
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
const ADB& b = rq_[ phase ].b;
|
||||
|
||||
@ -177,6 +177,8 @@ namespace Opm {
|
||||
computeRelPermWells(const SolutionState& state,
|
||||
const DataBlock& well_s,
|
||||
const std::vector<int>& well_cells) const;
|
||||
std::vector<ADB>
|
||||
computePressures(const SolutionState& state) const;
|
||||
|
||||
void
|
||||
computeMassFlux(const int actph ,
|
||||
|
||||
Loading…
Reference in New Issue
Block a user