From 13acc8ee03452725e4c4c00b6d47d86e9f22d2a8 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 19 Apr 2016 17:09:50 +0200 Subject: [PATCH] moving computeWellConnectionPressures to StandardWells the results look okay, while the running for flow_solvent needs further investigation even the results with flow_solvent actually look okay. With two different version of computePropertiesForWellConnectionPressures, flow_solvent produces the same results. This is something needs further investigation. The current implementation requires a copy of computeWellConnectionPressure in StandardWells and StandardWellsSolvent. That means probably we need to introduce the asImpl() for the Wells. --- opm/autodiff/BlackoilModelBase.hpp | 5 +- opm/autodiff/BlackoilModelBase_impl.hpp | 54 +++++-------------- opm/autodiff/BlackoilSolventModel.hpp | 2 +- opm/autodiff/StandardWells.hpp | 52 ++++++++++-------- opm/autodiff/StandardWellsSolvent.hpp | 11 ++++ opm/autodiff/StandardWellsSolvent_impl.hpp | 36 +++++++++++++ opm/autodiff/StandardWells_impl.hpp | 36 +++++++++++++ .../fullyimplicit/BlackoilPolymerModel.hpp | 2 +- .../BlackoilPolymerModel_impl.hpp | 4 +- 9 files changed, 135 insertions(+), 67 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8cf69a88..8d60bf619 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -380,9 +380,6 @@ namespace Opm { computeAccum(const SolutionState& state, const int aix ); - void computeWellConnectionPressures(const SolutionState& state, - const WellState& xw); - void assembleMassBalanceEq(const SolutionState& state); @@ -410,6 +407,8 @@ namespace Opm { const std::vector& mob_perfcells, const std::vector& b_perfcells, WellState& well_state); + + void addWellContributionToMassBalanceEq(const std::vector& cq_s, const SolutionState& state, const WellState& xw); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 5453d1d87..74cb2e980 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -798,41 +798,6 @@ namespace detail { - template - void - BlackoilModelBase:: - computeWellConnectionPressures(const SolutionState& state, - const WellState& xw) - { - if( ! localWellsActive() ) return ; - - using namespace Opm::AutoDiffGrid; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. - std::vector b_perf; - std::vector rsmax_perf; - std::vector rvmax_perf; - std::vector surf_dens_perf; - asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // Extract well connection depths. - const typename WellModel::Vector depth = cellCentroidsZToEigen(grid_); - const typename WellModel::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); - - // Gravity - const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - - asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav); - - } - - - - - template void BlackoilModelBase:: @@ -842,6 +807,9 @@ namespace detail { { using namespace Opm::AutoDiffGrid; + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = cellCentroidsZToEigen(grid_); + // If we have VFP tables, we need the well connection // pressures for the "simple" hydrostatic correction // between well depth and vfp table depth. @@ -849,14 +817,15 @@ namespace detail { SolutionState state = asImpl().variableState(reservoir_state, well_state); SolutionState state0 = state; asImpl().makeConstantState(state0); - asImpl().computeWellConnectionPressures(state0, well_state); + // asImpl().computeWellConnectionPressures(state0, well_state); + // Extract well connection depths. + asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells // asImpl().updateWellControls(well_state); // asImpl().stdWells().updateWellControls(well_state); - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); // Create the primary variables. @@ -869,7 +838,8 @@ namespace detail { // Compute initial accumulation contributions // and well connection pressures. asImpl().computeAccum(state0, 0); - asImpl().computeWellConnectionPressures(state0, well_state); + // asImpl().computeWellConnectionPressures(state0, well_state); + asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // OPM_AD_DISKVAL(state.pressure); @@ -1133,6 +1103,9 @@ namespace detail { } } + // gravity + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + int it = 0; bool converged; do { @@ -1147,7 +1120,6 @@ namespace detail { asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_); - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); converged = getWellConvergence(it); @@ -1201,7 +1173,9 @@ namespace detail { std::vector old_derivs = state.qs.derivative(); state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); } - asImpl().computeWellConnectionPressures(state, well_state); + // asImpl().computeWellConnectionPressures(state, well_state); + const ADB::V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); + asImpl().stdWells().computeWellConnectionPressures(state, well_state, fluid_, active_, phaseCondition(), depth, gravity); } if (!converged) { diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 1e0409667..d3b09322e 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -146,7 +146,7 @@ namespace Opm { using Base::drMaxRel; using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; + // using Base::computeWellConnectionPressures; // using Base::addWellControlEq; // using Base::computePropertiesForWellConnectionPressures; diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 0b1d7e591..c10f5ba85 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -133,29 +133,39 @@ namespace Opm { const VFPProperties& vfp_properties, WellState& well_state); - template - void updateWellControls(const Opm::PhaseUsage& pu, - const double gravity, - const VFPProperties& vfp_properties, - const bool terminal_output, - const std::vector& active, - WellState& xw) const; + template + void updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const; - // TODO: should LinearisedBlackoilResidual also be a template class? - template - void addWellFluxEq(const std::vector& cq_s, - const SolutionState& state, - LinearisedBlackoilResidual& residual); + // TODO: should LinearisedBlackoilResidual also be a template class? + template + void addWellFluxEq(const std::vector& cq_s, + const SolutionState& state, + LinearisedBlackoilResidual& residual); + + // TODO: some parameters, like gravity, maybe it is better to put in the member list + template + void addWellControlEq(const SolutionState& state, + const WellState& xw, + const Vector& aliveWells, + const std::vector active, + const VFPProperties& vfp_properties, + const double gravity, + LinearisedBlackoilResidual& residual); + + template + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity); - // TODO: some parameters, like gravity, maybe it is better to put in the member list - template - void addWellControlEq(const SolutionState& state, - const WellState& xw, - const Vector& aliveWells, - const std::vector active, - const VFPProperties& vfp_properties, - const double gravity, - LinearisedBlackoilResidual& residual); diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp index cc491d323..5801aadb0 100644 --- a/opm/autodiff/StandardWellsSolvent.hpp +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -34,6 +34,7 @@ namespace Opm { public: using Base = StandardWells; + using Base::computeWellConnectionDensitesPressures; // --------- Public methods --------- explicit StandardWellsSolvent(const Wells* wells); @@ -63,6 +64,16 @@ namespace Opm { const std::vector& active, std::vector& mob_perfcells, std::vector& b_perfcells) const; + + template + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity); + protected: const SolventPropsAdFromDeck* solvent_props_; int solvent_pos_; diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp index 3529e8708..fa5babbac 100644 --- a/opm/autodiff/StandardWellsSolvent_impl.hpp +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -197,6 +197,42 @@ namespace Opm + template + void + StandardWellsSolvent:: + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition, + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Extract well connection depths + // TODO: depth_perf should be a member of the StandardWells class + const Vector pdepth = subset(depth, wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + + + + template void StandardWellsSolvent:: diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 49efed601..2298d2661 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -289,6 +289,42 @@ namespace Opm + template + void + StandardWells:: + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition, + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Extract well connection depths + // TODO: depth_perf should be a member of the StandardWells class + const Vector pdepth = subset(depth, wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + + + + template void StandardWells:: diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 2f4c1a0e7..e63c80c58 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -202,7 +202,7 @@ namespace Opm { using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; + // using Base::computeWellConnectionPressures; // using Base::addWellControlEq; using Base::computeRelPerm; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index d6cb9b1cd..fce5c27b1 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -499,6 +499,7 @@ namespace Opm { // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = cellCentroidsZToEigen(grid_); // updateWellControls(well_state); stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); @@ -512,7 +513,8 @@ namespace Opm { // Compute initial accumulation contributions // and well connection pressures. computeAccum(state0, 0); - computeWellConnectionPressures(state0, well_state); + // computeWellConnectionPressures(state0, well_state); + stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // OPM_AD_DISKVAL(state.pressure);