From 977395fccd90c22ee1a007f5fb097c4e507aa9d9 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 26 Nov 2013 13:51:10 +0100 Subject: [PATCH] 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... --- opm/autodiff/BlackoilPropsAd.cpp | 59 +++++++++++++++++ opm/autodiff/BlackoilPropsAd.hpp | 13 ++++ opm/autodiff/BlackoilPropsAdFromDeck.cpp | 58 +++++++++++++++++ opm/autodiff/BlackoilPropsAdFromDeck.hpp | 13 ++++ opm/autodiff/BlackoilPropsAdInterface.hpp | 14 ++++ opm/autodiff/FullyImplicitBlackoilSolver.cpp | 67 ++++++++++++++++---- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 6 +- 7 files changed, 215 insertions(+), 15 deletions(-) diff --git a/opm/autodiff/BlackoilPropsAd.cpp b/opm/autodiff/BlackoilPropsAd.cpp index 0234f54b0..e41cb27e2 100644 --- a/opm/autodiff/BlackoilPropsAd.cpp +++ b/opm/autodiff/BlackoilPropsAd.cpp @@ -592,5 +592,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/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index fefd612df..88f5a30fd 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const; + private: const BlackoilPropertiesInterface& props_; PhaseUsage pu_; diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 9b4a57ad4..8902972cd 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -686,5 +686,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/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 64ad3a651..943a06f65 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const; + private: RockFromDeck rock_; boost::scoped_ptr satprops_; diff --git a/opm/autodiff/BlackoilPropsAdInterface.hpp b/opm/autodiff/BlackoilPropsAdInterface.hpp index ed798d717..ec98951bc 100644 --- a/opm/autodiff/BlackoilPropsAdInterface.hpp +++ b/opm/autodiff/BlackoilPropsAdInterface.hpp @@ -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 capPress(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const = 0; }; } // namespace Opm diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 260100bac..65e28396e 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -581,16 +581,17 @@ namespace { // for each active phase. const V transi = subset(geo_.transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - for (int phase = 0; phase < fluid_.numPhases(); ++phase) { - computeMassFlux(phase, transi, kr, state); + const std::vector 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 + FullyImplicitBlackoilSolver::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 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 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& 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 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 upwind(grid_, ops_, head.value()); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 815c8a4fd..042b2cb39 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -168,6 +168,9 @@ namespace Opm { BlackoilState& state, WellState& well_state) const; + std::vector + computePressures(const SolutionState& state) const; + std::vector computeRelPerm(const SolutionState& state) const; @@ -179,7 +182,8 @@ namespace Opm { void computeMassFlux(const int actph , const V& transi, - const std::vector& kr , + const ADB& kr , + const ADB& p , const SolutionState& state ); double