mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
include capillary pressure in the PDEs
I'm neither sure that this is fully correct nor that I found all occurences (so far, the output writing code is missing in this patch), but it seems to work for SPE1...
This commit is contained in:
Andreas Lauser
parent
96454a8496
commit
977395fccd
@ -592,5 +592,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
|
||||
|
||||
|
||||
@ -246,6 +246,19 @@ namespace Opm
|
||||
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_;
|
||||
PhaseUsage pu_;
|
||||
|
||||
@ -686,5 +686,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
|
||||
|
||||
|
||||
@ -247,6 +247,19 @@ namespace Opm
|
||||
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_;
|
||||
boost::scoped_ptr<SaturationPropsInterface> satprops_;
|
||||
|
||||
@ -251,6 +251,20 @@ 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
|
||||
|
||||
@ -581,16 +581,17 @@ namespace {
|
||||
// for each active phase.
|
||||
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
|
||||
computeMassFlux(phase, transi, kr, state);
|
||||
const std::vector<ADB> pressures = computePressures(state);
|
||||
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
|
||||
computeMassFlux(phaseIdx, transi, kr[phaseIdx], pressures[phaseIdx], state);
|
||||
// std::cout << "===== kr[" << phase << "] = \n" << std::endl;
|
||||
// std::cout << kr[phase];
|
||||
// std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
|
||||
// std::cout << rq_[phase].mflux;
|
||||
|
||||
residual_.mass_balance[ phase ] =
|
||||
pvdt*(rq_[phase].accum[1] - rq_[phase].accum[0])
|
||||
+ ops_.div*rq_[phase].mflux;
|
||||
residual_.mass_balance[ phaseIdx ] =
|
||||
pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
|
||||
+ ops_.div*rq_[phaseIdx].mflux;
|
||||
|
||||
// DUMP(residual_.mass_balance[phase]);
|
||||
}
|
||||
@ -968,6 +969,44 @@ namespace {
|
||||
}
|
||||
|
||||
|
||||
std::vector<ADB>
|
||||
FullyImplicitBlackoilSolver::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 Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const ADB sw = (active_[ Water ]
|
||||
? state.saturation[ pu.phase_pos[ Water ] ]
|
||||
: null);
|
||||
|
||||
const ADB so = (active_[ Oil ]
|
||||
? state.saturation[ pu.phase_pos[ Oil ] ]
|
||||
: null);
|
||||
|
||||
const ADB sg = (active_[ Gas ]
|
||||
? state.saturation[ pu.phase_pos[ Gas ] ]
|
||||
: null);
|
||||
|
||||
// convert the pressure offsets to the capillary pressures
|
||||
std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
|
||||
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
|
||||
#warning "what's the reference phase??"
|
||||
if (phaseIdx == BlackoilPhases::Liquid)
|
||||
continue;
|
||||
pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
|
||||
}
|
||||
|
||||
// add the total pressure to the capillary pressures
|
||||
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
|
||||
pressure[phaseIdx] += state.pressure;
|
||||
}
|
||||
|
||||
return pressure;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1005,30 +1044,30 @@ namespace {
|
||||
void
|
||||
FullyImplicitBlackoilSolver::computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const SolutionState& state )
|
||||
const ADB& kr ,
|
||||
const ADB& pressure,
|
||||
const SolutionState& state)
|
||||
{
|
||||
const int phase = canph_[ actph ];
|
||||
const int canonicalPhaseIdx = canph_[ actph ];
|
||||
|
||||
std::vector<PhasePresence> cond;
|
||||
classifyCondition(state, cond);
|
||||
|
||||
const ADB mu = fluidViscosity(phase, state.pressure, state.rs, cond, cells_);
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
const ADB mu = fluidViscosity(phase, pressure, state.rs, cond, cells_);
|
||||
|
||||
rq_[ actph ].mob = tr_mult * kr[ phase ] / mu;
|
||||
rq_[ actph ].mob = tr_mult * kr / mu;
|
||||
|
||||
const ADB rho = fluidDensity(phase, state.pressure, state.rs, cond, cells_);
|
||||
const ADB rho = fluidDensity(phase, pressure, state.rs, cond, cells_);
|
||||
|
||||
ADB& head = rq_[ actph ].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 * pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
|
||||
head = transi*dp;
|
||||
//head = transi*(ops_.ngrad * state.pressure) + gflux;
|
||||
//head = transi*(ops_.ngrad * pressure) + gflux;
|
||||
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
|
||||
|
||||
@ -168,6 +168,9 @@ namespace Opm {
|
||||
BlackoilState& state,
|
||||
WellState& well_state) const;
|
||||
|
||||
std::vector<ADB>
|
||||
computePressures(const SolutionState& state) const;
|
||||
|
||||
std::vector<ADB>
|
||||
computeRelPerm(const SolutionState& state) const;
|
||||
|
||||
@ -179,7 +182,8 @@ namespace Opm {
|
||||
void
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const ADB& kr ,
|
||||
const ADB& p ,
|
||||
const SolutionState& state );
|
||||
|
||||
double
|
||||
|
||||
Loading…
Reference in New Issue
Block a user