From 075e16dc360e72f4366beb288e24ef9a812bc020 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Mon, 27 Jan 2014 16:26:48 +0800 Subject: [PATCH] add capillary pressure for incom solver. --- ...FullyImplicitCompressiblePolymerSolver.cpp | 2 +- .../FullyImplicitTwophasePolymerSolver.cpp | 23 +++++++++- .../FullyImplicitTwophasePolymerSolver.hpp | 2 + .../fullyimplicit/IncompPropsAdFromDeck.cpp | 42 +++++++++++++++++++ .../fullyimplicit/IncompPropsAdFromDeck.hpp | 11 +++++ .../fullyimplicit/IncompPropsAdInterface.hpp | 12 ++++++ 6 files changed, 89 insertions(+), 3 deletions(-) diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 6566b7b2a..dfe6ddc09 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -831,7 +831,7 @@ namespace { // convert the pressure offsets to the capillary pressures std::vector pressure = fluid_.capPress(sw, so, sg, cells_); - pressure[0] = pressure[0] - pressure[0]; + pressure[0] = pressure[0] - pressure[1]; // add the total pressure to the capillary pressures for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) { diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index f24ac775a..5887bd947 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -548,6 +548,25 @@ namespace { } + std::vector + FullyImplicitTwophasePolymerSolver:: + computePressures(const SolutionState& state) const + { + const ADB sw = state.saturation[0]; + const ADB so = state.saturation[1]; + + // convert the pressure offsets to the capillary pressures + std::vector pressure = fluid_.capPress(sw, so, cells_); + pressure[0] = pressure[0] - pressure[1]; + + // add the total pressure to the capillary pressures + for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) { + pressure[phaseIdx] += state.pressure; + } + + return pressure; + } + void FullyImplicitTwophasePolymerSolver::computeMassFlux(const V& trans, const ADB& mc, @@ -561,18 +580,18 @@ namespace { rq_[1].mob = kro / V::Constant(kro.size(), 1, mus[1]); rq_[2].mob = mc * krw_eff * inv_wat_eff_vis; - const int nc = grid_.number_of_cells; V z(nc); // Compute z coordinates for (int c = 0; c < nc; ++c){ z[c] = grid_.cell_centroids[c * 3 + 2]; } + std::vector press = computePressures(state); for (int phase = 0; phase < 2; ++phase) { const ADB rho = fluidDensity(phase, state.pressure); ADB& head = rq_[phase].head; const ADB rhoavg = ops_.caver * rho; - const ADB dp = ops_.ngrad * state.pressure + const ADB dp = ops_.ngrad * press[phase] - gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix())); head = trans * dp; UpwindSelector upwind(grid_, ops_, head.value()); diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index 525b80a7d..9e0f5234a 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -102,6 +102,8 @@ namespace Opm { V transmissibility() const; + std::vector + computePressures(const SolutionState& state) const; void computeMassFlux(const V& trans, const ADB& mc, diff --git a/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp b/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp index 4781c85b8..4ca86247b 100644 --- a/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp +++ b/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp @@ -132,5 +132,47 @@ namespace Opm } return relperms; } + std::vector + IncompPropsAdFromDeck::capPress(const ADB& sw, + const ADB& so, + const Cells& cells) const + { + + const int n = cells.size(); + const int np = numPhases(); + const int num_blocks = so.numBlocks(); + Block s_all(n, np); + assert(sw.size() == n && so.size() == n); + s_all.col(0) = sw.value(); + s_all.col(1) = so.value(); + + Block pc(n, np); + Block dpc(n, np * np); + + satprops_.capPress(n, s_all.data(), cells.data(), pc.data(), dpc.data()); + + std::vector capPressures; + capPressures.reserve(2); + const ADB* s[2] = { &sw, &so}; + for (int phase1 = 0; phase1 < 3; ++phase1) { + const int phase1_pos = phase1; + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols()); + } + for (int phase2 = 0; phase2 < 3; ++phase2) { + const int phase2_pos = phase2; + // Assemble dpc1/ds2. + const int column = phase1_pos + phase2_pos; // Recall: Fortran ordering from props_.relperm() + ADB::M dpc1_ds2_diag = spdiag(dpc.col(column)); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block]; + } + } + capPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs)); + } + return capPressures; + } + } //namespace Opm diff --git a/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp b/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp index 8058611d0..c82a20685 100644 --- a/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp +++ b/opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp @@ -46,6 +46,17 @@ namespace Opm const ADB& so, 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] 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 Cells& cells) const; + private: RockFromDeck rock_; PvtPropertiesIncompFromDeck pvt_; diff --git a/opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp b/opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp index b74004c93..817cf7820 100644 --- a/opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp +++ b/opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp @@ -82,6 +82,18 @@ namespace Opm std::vector relperm(const ADB& sw, const ADB& so, 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] 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 Cells& cells) const = 0; + }; } // namespace Opm