From c4d567c5e499ed2743ec023808a4e2c2687439e1 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Mon, 27 Jan 2014 14:48:26 +0800 Subject: [PATCH] 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. --- opm/polymer/fullyimplicit/BlackoilPropsAd.cpp | 59 +++++++++++++++++++ opm/polymer/fullyimplicit/BlackoilPropsAd.hpp | 12 ++++ .../fullyimplicit/BlackoilPropsAdFromDeck.cpp | 58 ++++++++++++++++++ .../fullyimplicit/BlackoilPropsAdFromDeck.hpp | 13 ++++ .../BlackoilPropsAdInterface.hpp | 15 +++++ ...FullyImplicitCompressiblePolymerSolver.cpp | 31 +++++++++- ...FullyImplicitCompressiblePolymerSolver.hpp | 2 + 7 files changed, 188 insertions(+), 2 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPropsAd.cpp b/opm/polymer/fullyimplicit/BlackoilPropsAd.cpp index d83efc938..b09e21f7b 100644 --- a/opm/polymer/fullyimplicit/BlackoilPropsAd.cpp +++ b/opm/polymer/fullyimplicit/BlackoilPropsAd.cpp @@ -584,5 +584,64 @@ namespace Opm return relperms; } + std::vector 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 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 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 diff --git a/opm/polymer/fullyimplicit/BlackoilPropsAd.hpp b/opm/polymer/fullyimplicit/BlackoilPropsAd.hpp index b6fd925d5..f62793e84 100644 --- a/opm/polymer/fullyimplicit/BlackoilPropsAd.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPropsAd.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const; private: const BlackoilPropertiesInterface& props_; diff --git a/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp b/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp index c9ec833ab..4333a6a3c 100644 --- a/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp +++ b/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp @@ -678,5 +678,63 @@ namespace Opm return relperms; } + std::vector 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 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 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 diff --git a/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp b/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp index 5e0b983f4..d584ffb90 100644 --- a/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const; + private: RockFromDeck rock_; diff --git a/opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp b/opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp index a1417b7d6..153d993aa 100644 --- a/opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const = 0; + }; } // namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 24760f593..6566b7b2a 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -504,7 +504,6 @@ namespace { // for each active phase. const V trans = subset(geo_.transmissibility(), ops_.internal_faces); const std::vector 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 + FullyImplicitCompressiblePolymerSolver:: + computePressures(const SolutionState& state) const + { + const int nc = grid_.number_of_cells; + const std::vector& 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 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 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 upwind(grid_, ops_, head.value()); const ADB& b = rq_[ phase ].b; diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 3d2d44634..6d11d3da8 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -177,6 +177,8 @@ namespace Opm { computeRelPermWells(const SolutionState& state, const DataBlock& well_s, const std::vector& well_cells) const; + std::vector + computePressures(const SolutionState& state) const; void computeMassFlux(const int actph ,