From 5cd7468b51a621f9da49c5431de9c621f84d277b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 11 Aug 2016 13:31:52 +0200 Subject: [PATCH 1/8] Implement new well model Start using the well model desribed in SPE 12259 "Enhancements to the Strongly Coupled, Fully Implicit Well Model: Wellbore CrossFlow Modeling and Collective Well Control" The new well model uses three well primary variables: (the old one had 4) 1) bhp for rate controlled wells or total_rate for bhp controlled wells 2) flowing fraction of water Fw = g_w * Q_w / (g_w * Q_w + q_o * Q_o + q_g + Q_g) 3) flowing fraction of gas Fg = g_g * Q_g / (g_w * Q_w + q_o * Q_o + q_g + Q_g) where g_g = 0.01 and q_w = q_o = 1; Note 1: This is the starting point of implementing a well model using denseAD The plan is to gradually remove Eigen from the well model and instead use the fluidState object provided by the Ebos simulator Note 2: This is still work in progress and substantial cleaning and testing is still needed. Note 3: Runs SPE1, SPE9 and norne until 950 days. --- opm/autodiff/BlackoilModelEbos.hpp | 89 +- opm/autodiff/DefaultBlackoilSolutionState.hpp | 2 + .../NewtonIterationBlackoilInterleaved.cpp | 12 +- .../SimulatorFullyImplicitBlackoilEbos.hpp | 18 +- opm/autodiff/StandardWellsDense.hpp | 298 +++ opm/autodiff/StandardWellsDense_impl.hpp | 2237 +++++++++++++++++ .../WellStateFullyImplicitBlackoil.hpp | 82 + 7 files changed, 2684 insertions(+), 54 deletions(-) create mode 100644 opm/autodiff/StandardWellsDense.hpp create mode 100644 opm/autodiff/StandardWellsDense_impl.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index cc9100046..6ce28724a 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -28,7 +28,7 @@ #include #include -#include +#include #include #include #include @@ -41,6 +41,7 @@ #include #include #include +#include #include #include @@ -137,7 +138,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, - const StandardWells& well_model, + const StandardWellsDense& well_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output) : ebosSimulator_(ebosSimulator) @@ -300,6 +301,7 @@ namespace Opm { // Compute initial accumulation contributions // and well connection pressures. wellModel().computeWellConnectionPressures(state0, well_state); + wellModel().computeAccumWells(state0); } IterationReport iter_report = {false, false, 0, 0}; @@ -307,27 +309,27 @@ namespace Opm { return iter_report; } + double dt = timer.currentStepLength(); std::vector mob_perfcells; std::vector b_perfcells; wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); - if (param_.solve_welleq_initially_ && iterationIdx == 0) { + if (param_.solve_welleq_initially_ && iterationIdx == 0 ) { // solve the well equations as a pre-processing step - iter_report = solveWellEq(mob_perfcells, b_perfcells, state, well_state); + iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, state, well_state); } V aliveWells; std::vector cq_s; - wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + wellModel().computeWellFluxDense(state, mob_perfcells, b_perfcells, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); - wellModel().addWellFluxEq(cq_s, state, residual_); + wellModel().addWellFluxEq(cq_s, state, dt, residual_); addWellContributionToMassBalanceEq(cq_s, state, well_state); - wellModel().addWellControlEq(state, well_state, aliveWells, residual_); + //wellModel().addWellControlEq(state, well_state, aliveWells, residual_); if (param_.compute_well_potentials_) { SolutionState state0 = state; makeConstantState(state0); wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); } - return iter_report; } @@ -894,7 +896,7 @@ namespace Opm { std::vector phaseCondition_; // Well Model - StandardWells well_model_; + StandardWellsDense well_model_; V isRs_; V isRv_; @@ -915,8 +917,8 @@ namespace Opm { public: /// return the StandardWells object - StandardWells& wellModel() { return well_model_; } - const StandardWells& wellModel() const { return well_model_; } + StandardWellsDense& wellModel() { return well_model_; } + const StandardWellsDense& wellModel() const { return well_model_; } /// return the Well struct in the StandardWells const Wells& wells() const { return well_model_.wells(); } @@ -964,7 +966,7 @@ namespace Opm { std::vector indices = {{Pressure, Sw, Xvar}}; int foo = indices.size(); - indices.resize(5); + indices.resize(4); wellModel().variableStateWellIndices(indices, foo); const ADB& ones = ADB::constant(V::Ones(nc, 1)); @@ -1005,7 +1007,7 @@ namespace Opm { //state.saturation[Oil] = std::move(so); // wells - wellModel().variableStateExtractWellsVars(indices, vars, state); + wellModel().variableStateExtractWellsVars(indices, vars, state, xw); } V @@ -1218,14 +1220,14 @@ namespace Opm { for( int flowPhaseIdx = 0; flowPhaseIdx < numPhases; ++flowPhaseIdx ) { - adbJacs[ flowPhaseIdx ].resize( numPhases + 2 ); + adbJacs[ flowPhaseIdx ].resize( numPhases + 1 ); for( int pvIdx = 0; pvIdx < numPhases; ++pvIdx ) { jacs[ flowPhaseIdx ][ pvIdx ].finalize(); adbJacs[ flowPhaseIdx ][ pvIdx ].assign( std::move(jacs[ flowPhaseIdx ][ pvIdx ]) ); } // add two "dummy" matrices for the well primary variables - for( int pvIdx = numPhases; pvIdx < numPhases + 2; ++pvIdx ) { + for( int pvIdx = numPhases; pvIdx < numPhases + 1; ++pvIdx ) { adbJacs[ flowPhaseIdx ][ pvIdx ] = sparsityPattern.derivative()[pvIdx]; } @@ -1273,14 +1275,14 @@ namespace Opm { // create the Jacobian matrices for the legacy state. here we assume that the // sparsity pattern of the inputs is already correct /////// - std::vector poJac(numPhases + 2); + std::vector poJac(numPhases + 1); //std::vector TJac(numPhases + 2); std::vector> SJac(numPhases); std::vector> mobJac(numPhases); std::vector> bJac(numPhases); std::vector> pJac(numPhases); - std::vector RsJac(numPhases + 2); - std::vector RvJac(numPhases + 2); + std::vector RsJac(numPhases + 1); + std::vector RvJac(numPhases + 1); // reservoir stuff for (int pvIdx = 0; pvIdx < numPhases; ++ pvIdx) { @@ -1296,7 +1298,7 @@ namespace Opm { } // auxiliary equations - for (int pvIdx = numPhases; pvIdx < numPhases + 2; ++ pvIdx) { + for (int pvIdx = numPhases; pvIdx < numPhases + 1; ++ pvIdx) { legacyState.pressure.derivative()[pvIdx].toSparse(poJac[pvIdx]); //legacyState.temperature.derivative()[pvIdx].toSparse(TJac[pvIdx]); legacyState.rs.derivative()[pvIdx].toSparse(RsJac[pvIdx]); @@ -1304,10 +1306,10 @@ namespace Opm { } for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - SJac[phaseIdx].resize(numPhases + 2); - mobJac[phaseIdx].resize(numPhases + 2); - bJac[phaseIdx].resize(numPhases + 2); - pJac[phaseIdx].resize(numPhases + 2); + SJac[phaseIdx].resize(numPhases + 1); + mobJac[phaseIdx].resize(numPhases + 1); + bJac[phaseIdx].resize(numPhases + 1); + pJac[phaseIdx].resize(numPhases + 1); for (int pvIdx = 0; pvIdx < numPhases; ++ pvIdx) { SJac[phaseIdx][pvIdx].resize(numCells, numCells); SJac[phaseIdx][pvIdx].reserve(numCells); @@ -1323,7 +1325,7 @@ namespace Opm { } // auxiliary equations for the saturations and pressures - for (int pvIdx = numPhases; pvIdx < numPhases + 2; ++ pvIdx) { + for (int pvIdx = numPhases; pvIdx < numPhases + 1; ++ pvIdx) { legacyState.saturation[phaseIdx].derivative()[pvIdx].toSparse(SJac[phaseIdx][pvIdx]); legacyState.saturation[phaseIdx].derivative()[pvIdx].toSparse(mobJac[phaseIdx][pvIdx]); legacyState.saturation[phaseIdx].derivative()[pvIdx].toSparse(bJac[phaseIdx][pvIdx]); @@ -1403,10 +1405,10 @@ namespace Opm { std::vector RsAdbJacs; std::vector RvAdbJacs; - poAdbJacs.resize(numPhases + 2); - RsAdbJacs.resize(numPhases + 2); - RvAdbJacs.resize(numPhases + 2); - for(int pvIdx = 0; pvIdx < numPhases + 2; ++pvIdx) + poAdbJacs.resize(numPhases + 1); + RsAdbJacs.resize(numPhases + 1); + RvAdbJacs.resize(numPhases + 1); + for(int pvIdx = 0; pvIdx < numPhases + 1; ++pvIdx) { poAdbJacs[pvIdx].assign(poJac[pvIdx]); RsAdbJacs[pvIdx].assign(RsJac[pvIdx]); @@ -1419,11 +1421,11 @@ namespace Opm { std::vector> pAdbJacs(numPhases); for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - SAdbJacs[phaseIdx].resize(numPhases + 2); - mobAdbJacs[phaseIdx].resize(numPhases + 2); - bAdbJacs[phaseIdx].resize(numPhases + 2); - pAdbJacs[phaseIdx].resize(numPhases + 2); - for(int pvIdx = 0; pvIdx < numPhases + 2; ++pvIdx) + SAdbJacs[phaseIdx].resize(numPhases + 1); + mobAdbJacs[phaseIdx].resize(numPhases + 1); + bAdbJacs[phaseIdx].resize(numPhases + 1); + pAdbJacs[phaseIdx].resize(numPhases + 1); + for(int pvIdx = 0; pvIdx < numPhases + 1; ++pvIdx) { SAdbJacs[phaseIdx][pvIdx].assign(SJac[phaseIdx][pvIdx]); mobAdbJacs[phaseIdx][pvIdx].assign(mobJac[phaseIdx][pvIdx]); @@ -1529,6 +1531,7 @@ namespace Opm { IterationReport solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, + const double dt, SolutionState& state, WellState& well_state) { @@ -1557,16 +1560,16 @@ namespace Opm { do { // bhp and Q for the wells std::vector vars0; - vars0.reserve(2); + vars0.reserve(1); wellModel().variableWellStateInitials(well_state, vars0); std::vector vars = ADB::variables(vars0); SolutionState wellSolutionState = state0; - wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState); - wellModel().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState, well_state); + wellModel().computeWellFluxDense(wellSolutionState, mob_perfcells_const, b_perfcells_const, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); - wellModel().addWellFluxEq(cq_s, wellSolutionState, residual_); - wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); + wellModel().addWellFluxEq(cq_s, wellSolutionState, dt, residual_); + //wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); converged = getWellConvergence(it); if (converged) { @@ -1577,9 +1580,9 @@ namespace Opm { if( localWellsActive() ) { std::vector eqs; - eqs.reserve(2); + eqs.reserve(1); eqs.push_back(residual_.well_flux_eq); - eqs.push_back(residual_.well_eq); + //eqs.push_back(residual_.well_eq); ADB total_residual = vertcatCollapseJacs(eqs); const std::vector& Jn = total_residual.derivative(); typedef Eigen::SparseMatrix Sp; @@ -1615,7 +1618,11 @@ namespace Opm { std::vector old_derivs = state.qs.derivative(); state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); } - computeWellConnectionPressures(state, well_state); + + const ADB::V new_well_var = Eigen::Map(& well_state.wellSolutions()[0], well_state.wellSolutions().size()); + std::vector old_derivs = state.wellVariables.derivative(); + state.wellVariables = ADB::function(std::move(new_well_var), std::move(old_derivs)); + //computeWellConnectionPressures(state, well_state); } if (!converged) { diff --git a/opm/autodiff/DefaultBlackoilSolutionState.hpp b/opm/autodiff/DefaultBlackoilSolutionState.hpp index 6ec69746d..06d7a4aef 100644 --- a/opm/autodiff/DefaultBlackoilSolutionState.hpp +++ b/opm/autodiff/DefaultBlackoilSolutionState.hpp @@ -37,6 +37,7 @@ namespace Opm { , rv ( ADB::null()) , qs ( ADB::null()) , bhp ( ADB::null()) + , wellVariables ( ADB::null()) , canonical_phase_pressures(3, ADB::null()) { } @@ -47,6 +48,7 @@ namespace Opm { ADB rv; ADB qs; ADB bhp; + ADB wellVariables; // Below are quantities stored in the state for optimization purposes. std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active. }; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 9baba698f..1d96d6ea3 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -390,7 +390,7 @@ namespace Opm //const int np = residual.material_balance_eq.size(); assert( np == int(residual.material_balance_eq.size()) ); std::vector eqs; - eqs.reserve(np + 2); + eqs.reserve(np + 1); for (int phase = 0; phase < np; ++phase) { eqs.push_back(residual.material_balance_eq[phase]); } @@ -401,14 +401,14 @@ namespace Opm if( hasWells ) { eqs.push_back(residual.well_flux_eq); - eqs.push_back(residual.well_eq); + //eqs.push_back(residual.well_eq); // Eliminate the well-related unknowns, and corresponding equations. - elim_eqs.reserve(2); + elim_eqs.reserve(1); elim_eqs.push_back(eqs[np]); eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns. - elim_eqs.push_back(eqs[np]); - eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. + //elim_eqs.push_back(eqs[np]); + //eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. assert(int(eqs.size()) == np); } @@ -500,7 +500,7 @@ namespace Opm if ( hasWells ) { // Compute full solution using the eliminated equations. // Recovery in inverse order of elimination. - dx = recoverVariable(elim_eqs[1], dx, np); + //dx = recoverVariable(elim_eqs[1], dx, np); dx = recoverVariable(elim_eqs[0], dx, np); } return dx; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 048b11b81..e52b573e3 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -25,11 +25,15 @@ #include #include #include +#include +#include + + namespace Opm { class SimulatorFullyImplicitBlackoilEbos; -class StandardWells; +class StandardWellsDense; /// a simulator for the blackoil model class SimulatorFullyImplicitBlackoilEbos @@ -44,7 +48,7 @@ public: typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; - typedef StandardWells WellModel; + typedef StandardWellsDense WellModel; /// Initialise from parameters and objects to observe. @@ -169,11 +173,11 @@ public: // -1 means that we'll take the last report step that was written const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); - output_writer_.restore( timer, - state, - prev_well_state, - restorefilename, - desiredRestoreStep ); +// output_writer_.restore( timer, +// state, +// prev_well_state, +// restorefilename, +// desiredRestoreStep ); } unsigned int totalLinearizations = 0; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp new file mode 100644 index 000000000..601760fdb --- /dev/null +++ b/opm/autodiff/StandardWellsDense.hpp @@ -0,0 +1,298 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_STANDARDWELLSDENSE_HEADER_INCLUDED +#define OPM_STANDARDWELLSDENSE_HEADER_INCLUDED + +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm { + + + /// Class for handling the standard well model. + class StandardWellsDense { + public: + struct WellOps { + explicit WellOps(const Wells* wells); + Eigen::SparseMatrix w2p; // well -> perf (scatter) + Eigen::SparseMatrix p2w; // perf -> well (gather) + std::vector well_cells; // the set of perforated cells + }; + + // --------- Types --------- + using ADB = AutoDiffBlock; + + typedef DenseAd::Evaluation Eval; + //typedef AutoDiffBlock ADB; + using Vector = ADB::V; + using V = ADB::V; + + // copied from BlackoilModelBase + // should put to somewhere better + using DataBlock = Eigen::Array; + // --------- Public methods --------- + explicit StandardWellsDense(const Wells* wells_arg); + + void init(const BlackoilPropsAdInterface* fluid_arg, + const std::vector* active_arg, + const std::vector* pc_arg, + const VFPProperties* vfp_properties_arg, + const double gravity_arg, + const Vector& depth_arg); + + const WellOps& wellOps() const; + + int numPhases() const { return wells().number_of_phases; }; + + const Wells& wells() const; + + const Wells* wellsPointer() const; + + /// return true if wells are available in the reservoir + bool wellsActive() const; + void setWellsActive(const bool wells_active); + /// return true if wells are available on this process + bool localWellsActive() const; + + int numWellVars() const; + + /// Density of each well perforation + Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel + const Vector& wellPerforationDensities() const; + + /// Diff to bhp for each well perforation. + Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel + const Vector& wellPerforationPressureDiffs() const; + + template + void extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + std::vector& mob_perfcells, + std::vector& b_perfcells) const; + + template + void computeWellFlux(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const; + + template + void + computeWellFluxDense(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + std::vector& cq_s) const; + + Eval extractDenseAD(const ADB& data, int i, int j) const; + Eval extractDenseADWell(const ADB& data, int i) const; + + const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; + + + template + void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const; + + template + void updateWellState(const Vector& dwells, + const double dpmaxrel, + WellState& well_state); + + template + void updateWellControls(WellState& xw) const; + + // TODO: should LinearisedBlackoilResidual also be a template class? + template + void addWellFluxEq(const std::vector& cq_s, + const SolutionState& state, + const double dt, + 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, + LinearisedBlackoilResidual& residual); + + template + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw); + + // state0 is non-constant, while it will not be used outside of the function + template + void + computeWellPotentials(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state0, + WellState& well_state); + + template + void + variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state, + WellState& well_state) const; + + + void + variableStateWellIndices(std::vector& indices, + int& next) const; + + std::vector + variableWellStateIndices() const; + + template + void + variableWellStateInitials(const WellState& xw, + std::vector& vars0) const; + + /// If set, computeWellFlux() will additionally store the + /// total reservoir volume perforation fluxes. + void setStoreWellPerforationFluxesFlag(const bool store_fluxes); + + /// Retrieves the stored fluxes. It is an error to call this + /// unless setStoreWellPerforationFluxesFlag(true) has been + /// called. + const Vector& getStoredWellPerforationFluxes() const; + + /// upate the dynamic lists related to economic limits + template + void + updateListEconLimited(ScheduleConstPtr schedule, + const int current_step, + const Wells* wells, + const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const; + + template + std::vector wellVolumeFractions(const SolutionState& state) const; + + + template + void computeAccumWells(const SolutionState& state); + + protected: + bool wells_active_; + const Wells* wells_; + const WellOps wops_; + + const BlackoilPropsAdInterface* fluid_; + const std::vector* active_; + const std::vector* phase_condition_; + const VFPProperties* vfp_properties_; + double gravity_; + // the depth of the all the cell centers + // for standard Wells, it the same with the perforation depth + Vector perf_cell_depth_; + + Vector well_perforation_densities_; + Vector well_perforation_pressure_diffs_; + + bool store_well_perforation_fluxes_; + Vector well_perforation_fluxes_; + + std::vector F0_; + + // protected methods + template + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav); + + + template + bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const int well_number) const; + + using WellMapType = typename WellState::WellMapType; + using WellMapEntryType = typename WellState::mapentry_t; + + // a tuple type for ratio limit check. + // first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three + // values should not be used. + // second value indicates whehter there is only one connection left. + // third value indicates the indx of the worst-offending connection. + // the last value indicates the extent of the violation for the worst-offending connection, which is defined by + // the ratio of the actual value to the value of the violated limit. + using RatioCheckTuple = std::tuple; + + enum ConnectionIndex { + INVALIDCONNECTION = -10000 + }; + + + template + RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const; + + template + RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const; + + }; + + +} // namespace Opm + +#include "StandardWellsDense_impl.hpp" + +#endif diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp new file mode 100644 index 000000000..46fd6b85d --- /dev/null +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -0,0 +1,2237 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + Copyright 2016 IRIS AS. + + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#include +#include + +#include +#include +#include +#include + + + + + +namespace Opm +{ + + + StandardWellsDense:: + WellOps::WellOps(const Wells* wells) + : w2p(), + p2w(), + well_cells() + { + if( wells ) + { + w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); + + const int nw = wells->number_of_wells; + const int* const wpos = wells->well_connpos; + + typedef Eigen::Triplet Tri; + + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } + } + + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + + well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); + } + } + + + + + + StandardWellsDense::StandardWellsDense(const Wells* wells_arg) + : wells_active_(wells_arg!=nullptr) + , wells_(wells_arg) + , wops_(wells_arg) + , fluid_(nullptr) + , active_(nullptr) + , phase_condition_(nullptr) + , vfp_properties_(nullptr) + , well_perforation_densities_(Vector()) + , well_perforation_pressure_diffs_(Vector()) + , store_well_perforation_fluxes_(false) + { + } + + + + + + void + StandardWellsDense::init(const BlackoilPropsAdInterface* fluid_arg, + const std::vector* active_arg, + const std::vector* pc_arg, + const VFPProperties* vfp_properties_arg, + const double gravity_arg, + const Vector& depth_arg) + { + fluid_ = fluid_arg; + active_ = active_arg; + phase_condition_ = pc_arg; + vfp_properties_ = vfp_properties_arg; + gravity_ = gravity_arg; + perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);; + } + + + + + + const Wells& StandardWellsDense::wells() const + { + assert(wells_ != 0); + return *(wells_); + } + + + const Wells* StandardWellsDense::wellsPointer() const + { + return wells_; + } + + + + bool StandardWellsDense::wellsActive() const + { + return wells_active_; + } + + + + + + void StandardWellsDense::setWellsActive(const bool wells_active) + { + wells_active_ = wells_active; + } + + + + + + bool StandardWellsDense::localWellsActive() const + { + return wells_ ? (wells_->number_of_wells > 0 ) : false; + } + + + + + + int + StandardWellsDense::numWellVars() const + { + if ( !localWellsActive() ) + { + return 0; + } + + // For each well, we have a bhp variable, and one flux per phase. + const int nw = wells().number_of_wells; + return (numPhases() + 1) * nw; + } + + + + + + const StandardWellsDense::WellOps& + StandardWellsDense::wellOps() const + { + return wops_; + } + + + + + + StandardWellsDense::Vector& StandardWellsDense::wellPerforationDensities() + { + return well_perforation_densities_; + } + + + + + + const StandardWellsDense::Vector& + StandardWellsDense::wellPerforationDensities() const + { + return well_perforation_densities_; + } + + + + + + StandardWellsDense::Vector& + StandardWellsDense::wellPerforationPressureDiffs() + { + return well_perforation_pressure_diffs_; + } + + + + + + const StandardWellsDense::Vector& + StandardWellsDense::wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } + + + + + template + void + StandardWellsDense:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + // const std::vector& pc = phaseCondition(); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = (*phase_condition_)[well_cells[perf]]; + } + const PhaseUsage& pu = fluid_->phaseUsage(); + DataBlock b(nperf, pu.num_phases); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + assert((*active_)[Oil]); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null(); + const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null(); + const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + } + if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { + const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + + const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } + + // b is row major, so can just copy data. + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + + + + + template + void + StandardWellsDense:: + computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav) + { + // Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid_->phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const int nperf = wells().well_connpos[wells().number_of_wells]; + + // Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, cd, grav); + + // Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf); + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + } + + + + + + template + void + StandardWellsDense:: + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + 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, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const Vector& pdepth = perf_cell_depth_; + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_); + + } + + + + + + template + void + StandardWellsDense:: + extractWellPerfProperties(const SolutionState& /* state */, + const std::vector& rq, + std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + // If we have wells, extract the mobilities and b-factors for + // the well-perforated cells. + if ( !localWellsActive() ) { + mob_perfcells.clear(); + b_perfcells.clear(); + return; + } else { + const std::vector& well_cells = wellOps().well_cells; + const int num_phases = wells().number_of_phases; + mob_perfcells.resize(num_phases, ADB::null()); + b_perfcells.resize(num_phases, ADB::null()); + for (int phase = 0; phase < num_phases; ++phase) { + mob_perfcells[phase] = subset(rq[phase].mob, well_cells); + b_perfcells[phase] = subset(rq[phase].b, well_cells); + } + } + } + + template + void + StandardWellsDense:: + computeWellFluxDense(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + Vector Tw = Eigen::Map(wells().WI, nperf); + const std::vector& well_cells = wellOps().well_cells; + std::vector well_id(nperf); + + // pressure diffs computed already (once per step, not changing per iteration) + const Vector& cdp = wellPerforationPressureDiffs(); + + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + std::vector cmix_s_ADB = wellVolumeFractions(state); + + const int oilpos = pu.phase_pos[Oil]; + ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; + const ADB rsSat = fluid_->rsSat(perfpressure, cmix_s_ADB[oilpos], well_cells); + const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s_ADB[oilpos], well_cells); + + for (int w = 0; w < nw; ++w) { + Eval bhp = extractDenseADWell(state.bhp,w); + + + // TODO: fix for 2-phase case + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + cmix_s[phase] = extractDenseADWell(cmix_s_ADB[phase],w); + } + + //std::cout <<"cmix gas "<< w<< " "< b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + b_perfcells_dense[phase] = extractDenseAD(b_perfcells[phase], perf, cell_idx); + mob_perfcells_dense[phase] = extractDenseAD(mob_perfcells[phase], perf, cell_idx); + + } + + // Pressure drawdown (also used to determine direction of flow) + Eval drawdown = pressure - bhp - cdp[perf]; + + // injection perforations + if ( drawdown.value > 0 ) { + + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) + continue; + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + const Eval cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + } + + + if ((*active_)[Oil] && (*active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const Eval cq_psOil = cq_ps[oilpos]; + const Eval cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs * cq_psOil; + cq_ps[oilpos] += rv * cq_psGas; + } + + + // map to ADB + for (int phase = 0; phase < np; ++phase) { + cq_s_dense[phase][perf] = cq_ps[phase]; + } + + } else { + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) + continue; + + // Using total mobilities + Eval total_mob_dense = mob_perfcells_dense[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob_dense += mob_perfcells_dense[phase]; + } + // injection perforations total volume rates + const Eval cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + Eval volumeRatio = 0.0; + if ((*active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if ((*active_)[Oil] && (*active_)[Gas]) { + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + Eval rvPerf = 0.0; + if (cmix_s[gaspos] > 0) + rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; + + if (rvPerf.value > rvSat.value()[w]) { + rvPerf = 0.0; + rvPerf.value = rvSat.value()[w]; + } + + Eval rsPerf = 0.0; + if (cmix_s[oilpos] > 0) + rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; + + if (rsPerf.value > rsSat.value()[w]) { + rsPerf = 0.0; + rsPerf.value = rsSat.value()[w]; + } + + // Incorporate RS/RV factors if both oil and gas active + const Eval d = 1.0 - rvPerf * rsPerf; + + const Eval tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < cq_s2; + //Vector aliveWells; + //computeWellFlux(state,mob_perfcells,b_perfcells, aliveWells,cq_s2); + + //for (int phase = 0; phase < np; ++phase) { + //if( !(((cq_s[phase].value() - cq_s2[phase].value()).abs()<1e-10).all()) ) { + //std::cout << "phase " << phase << std::endl; + //std::cout << cq_s2[phase].value() << std::endl; + //std::cout << cq_s[phase].value() << std::endl; + //} + //} + } + + + + + + + template + void + StandardWellsDense:: + computeWellFlux(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector Tw = Eigen::Map(wells().WI, nperf); + const std::vector& well_cells = wellOps().well_cells; + + // pressure diffs computed already (once per step, not changing per iteration) + const Vector& cdp = wellPerforationPressureDiffs(); + // Extract needed quantities for the perforation cells + const ADB& p_perfcells = subset(state.pressure, well_cells); + + // Perforation pressure + const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; + + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcells - perfpressure; + + // Compute vectors with zero and ones that + // selects the wanted quantities. + + // selects injection perforations + Vector selectInjectingPerforations = Vector::Zero(nperf); + // selects producing perforations + Vector selectProducingPerforations = Vector::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (drawdown.value()[c] < 0) + selectInjectingPerforations[c] = 1; + else + selectProducingPerforations[c] = 1; + } + std::vector distri(np); + for (int p = 0; p < np; ++p) { + distri[p] = V::Zero(nw); + } + + + Vector isResv = Vector::Zero(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const double* distr = well_controls_get_current_distr(wc); + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + isResv[w] = 1; + } + for (int p = 0; p < np; ++p) { + distri[p][w] = distr[p]; + //std::cout << "distr " << distr[p] << " comp_frac " << comp_frac << std::endl; + } + + } + + + // Handle cross flow + const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); + const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); + for (int w = 0; w < nw; ++w) { + if (!wells().allow_cf[w]) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // Crossflow is not allowed; reverse flow is prevented. + // At least one of the perforation must be open in order to have a meeningful + // equation to solve. For the special case where all perforations have reverse flow, + // and the target rate is non-zero all of the perforations are keept open. + if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { + selectProducingPerforations[perf] = 0.0; + } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ + selectInjectingPerforations[perf] = 0.0; + } + } + } + } + + // HANDLE FLOW INTO WELLBORE + // compute phase volumetric rates at standard conditions + std::vector cq_p(np, ADB::null()); + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); + cq_ps[phase] = b_perfcells[phase] * cq_p[phase]; + } + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + if ((*active_)[Oil] && (*active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const ADB cq_psOil = cq_ps[oilpos]; + const ADB cq_psGas = cq_ps[gaspos]; + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); + cq_ps[gaspos] += rs_perfcells * cq_psOil; + cq_ps[oilpos] += rv_perfcells * cq_psGas; + } + + // HANDLE FLOW OUT FROM WELLBORE + // Using total mobilities + ADB total_mob = mob_perfcells[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob += mob_perfcells[phase]; + } + // injection perforations total volume rates + const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); + + // Store well perforation total fluxes (reservor volumes) if requested. + if (store_well_perforation_fluxes_) { + // Ugly const-cast, but unappealing alternatives. + Vector& wf = const_cast(well_perforation_fluxes_); + wf = cqt_i.value(); + for (int phase = 0; phase < np; ++phase) { + wf += cq_p[phase].value(); + } + } + + // compute wellbore mixture for injecting perforations + // The wellbore mixture depends on the inflow from the reservoar + // and the well injection rates. + + // compute avg. and total wellbore phase volumetric rates at standard conds +// const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); +// std::vector wbq(np, ADB::null()); +// ADB wbqt = ADB::constant(Vector::Zero(nw)); +// for (int phase = 0; phase < np; ++phase) { +// const ADB& q_ps = wellOps().p2w * cq_ps[phase]; +// const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); +// Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); +// const int pos = pu.phase_pos[phase]; +// wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps; +// wbqt += wbq[phase]; +// } + // compute wellbore mixture at standard conditions. + + + std::vector cmix_s(np, ADB::null()); + Vector ones = Vector::Constant(nperf,1.0); + cmix_s[Water] = wellOps().w2p * subset(state.wellVariables, Span(nw, 1, 1*nw)); + cmix_s[Gas] = wellOps().w2p * subset(state.wellVariables, Span(nw, 1, 2*nw)); + cmix_s[Oil] = (ones - cmix_s[Water] - cmix_s[Gas]); +// std::cout << cmix_s[Gas].value() << std::endl; + +// std::vector g(np, Vector::Constant(nw,1.0)); +// g[Gas] = Vector::Constant(nw,0.01); +// const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); +// std::vector wbq(np, ADB::null()); +// ADB wbqt = ADB::constant(Vector::Zero(nw)); +// for (int phase = 0; phase < np; ++phase) { +// const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); +// wbq[phase] = q_s * ( (Vector::Constant(nw,1.0) - isResv) * g[phase] + isResv * distri[phase]); +// wbqt += wbq[phase]; +// } +// Selector notDeadWells_selector(wbqt.value(), Selector::Zero); +// std::vector cmix_s(np, ADB::null()); +// for (int phase = 0; phase < np; ++phase) { +// const int pos = pu.phase_pos[phase]; +// cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); +// } + + + // compute volume ratio between connection at standard conditions + ADB volumeRatio = ADB::constant(Vector::Zero(nperf)); + + if ((*active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells[watpos]; + } + + + if ((*active_)[Oil] && (*active_)[Gas]) { + // Incorporate RS/RV factors if both oil and gas active +// const ADB& rv = subset(state.rv, well_cells); +// const ADB& rs = subset(state.rs, well_cells); +// const ADB d = Vector::Constant(nperf,1.0) - rv* rs; + + + + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + + Selector noGas_selector(cmix_s[gaspos].value(), Selector::GreaterZero); + Selector noOil_selector(cmix_s[oilpos].value(), Selector::GreaterZero); + ADB rv = noGas_selector.select(cmix_s[oilpos]/cmix_s[gaspos], ADB::constant(Vector::Constant(nperf,0.0))); + ADB rs = noOil_selector.select(cmix_s[gaspos]/cmix_s[oilpos], ADB::constant(Vector::Constant(nperf,0.0))); + + const ADB rsSat = fluid_->rsSat(perfpressure, cmix_s[oilpos], well_cells); + const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s[oilpos], well_cells); + + Selector maxRs_selector(rs.value() - rsSat.value(), Selector::GreaterZero); + Selector maxRv_selector(rv.value() - rvSat.value(), Selector::GreaterZero); + rv = maxRv_selector.select(rvSat, rv); + rs = maxRs_selector.select(rsSat, rs); + + const ADB d = Vector::Constant(nperf,1.0) - rv * rs; + + + const ADB tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; + volumeRatio += tmp_oil / b_perfcells[oilpos]; + + const ADB tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; + volumeRatio += tmp_gas / b_perfcells[gaspos]; + } + else { + if ((*active_)[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos]; + } + if ((*active_)[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos]; + } + } + + + // injecting connections total volumerates at standard conditions + //std::cout << volumeRatio.value() << std::endl; + ADB cqt_is = cqt_i/volumeRatio; + + // connection phase volumerates at standard conditions + cq_s.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; //*b_perfcells[phase]; + } + + // check for dead wells (used in the well controll equations) + aliveWells = Vector::Constant(nw, 1.0); +// for (int w = 0; w < nw; ++w) { +// if (wbqt.value()[w] == 0) { +// aliveWells[w] = 0.0; +// } +// } + } + + typedef DenseAd::Evaluation Eval; + Eval + StandardWellsDense:: + extractDenseAD(const ADB& data, int i, int j) const + { + Eval output = 0.0; + output.value = data.value()[i]; + const int np = wells().number_of_phases; + const std::vector& jac = data.derivative(); + //std::cout << jac.size() << std::endl; + int numblocs = jac.size(); + for (int b = 0; b < numblocs; ++b) { + if (b < np) { // don't copy well blocks) + //std::cout << jac[b].coeff(i,j) << std::endl; + output.derivatives[b] = jac[b].coeff(i,j); + } + } + return output; + } + + typedef DenseAd::Evaluation Eval; + Eval + StandardWellsDense:: + extractDenseADWell(const ADB& data, int i) const + { + Eval output = 0.0; + output.value = data.value()[i]; + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const std::vector& jac = data.derivative(); + //std::cout << jac.size() << std::endl; + int numblocs = jac.size(); + for (int b = 0; b < np; ++b) { + output.derivatives[b+np] = jac[numblocs-1].coeff(i, b*nw + i); + } + return output; + } + + const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const + { + typedef typename ADB::M M; + const int nLocal = local.size(); + typename ADB::V value( nLocal ); + //const int numVars = 5; + const int np = wells().number_of_phases; + + std::vector> mat(np, Eigen::SparseMatrix(nLocal,nc)); + Eigen::SparseMatrix matFlux(nLocal,np*nw); + Eigen::SparseMatrix matBHP(nLocal,nw); + + for( int i=0; i jacs( numVars ); + if (numVars == 4) { + for( int d=0; d(deri[d]); + jacs[ d ] = M(mat[d]); + } + + jacs[3] = M(matFlux); + //jacs[4] = M(matBHP); + } + else if (numVars == 1) { + jacs[0] = M(matFlux); + //jacs[1] = M(matBHP); + } + //std::cout << numVars << std::endl; + + return ADB::function( std::move( value ), std::move( jacs )); + } + + + + + + template + void + StandardWellsDense:: + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); + + // Update the perforation pressures. + const Vector& cdp = wellPerforationPressureDiffs(); + const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp; + xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); + } + + + + + + + + template + void + StandardWellsDense:: + updateWellState(const Vector& dwells, + const double dpmaxrel, + WellState& well_state) + { + if( localWellsActive() ) + { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const Vector dxvar_well = subset(dwells, Span(np*nw, 1, varstart)); + //const Vector dqs = subset(dwells, Span(np*nw, 1, varstart)); + varstart += dxvar_well.size(); + //const Vector dbhp = subset(dwells, Span(nw, 1, varstart)); + //varstart += dbhp.size(); + assert(varstart == dwells.size()); + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const Vector xvar_well_old = Eigen::Map(&well_state.wellSolutions()[0], nw*np); + + double dFLimit = 0.2; + double dBHPLimit = 2; + double dTotalRateLimit = 0.5; + //std::cout << "dxvar_well "< F(np,0.0); + const int sign2 = dxvar_well[nw + w] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dxvar_well[nw + w]),dFLimit); + well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited; + const int sign3 = dxvar_well[2*nw + w] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dxvar_well[2*nw + w]),dFLimit); + well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited; + F[Water] = well_state.wellSolutions()[nw + w]; + F[Gas] = well_state.wellSolutions()[2*nw + w]; + F[Oil] = 1.0 - F[Water] - F[Gas]; + +// const double dFw = dxvar_well[nw + w]; +// const double dFg = dxvar_well[nw*2 + w]; +// const double dFo = - dFw - dFg; +// //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; +// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; +// step = std::min(step, 1.0); +// //std::cout << step << std::endl; +// F[Water] = xvar_well_old[nw + w] - step*dFw; +// F[Gas] = xvar_well_old[2*nw + w] - step*dFg; +// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; + + if (F[Water] < 0.0) { + F[Gas] /= (1.0 - F[Water]); + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + if (F[Gas] < 0.0) { + F[Water] /= (1.0 - F[Gas]); + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } + if (F[Oil] < 0.0) { + F[Water] /= (1.0 - F[Oil]); + F[Gas] /= (1.0 - F[Oil]); + F[Oil] = 0.0; + } + well_state.wellSolutions()[nw + w] = F[Water]; + well_state.wellSolutions()[2*nw + w] = F[Gas]; + + //std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + + std::vector g = {1,1,0.01}; + + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + for (int p = 0; p < np; ++p) { + F[p] /= distr[p]; + } + } else { + for (int p = 0; p < np; ++p) { + F[p] /= g[p]; + } + } + //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + + + + + +// const double dFw = dxvar_well[nw + w]; +// const double dFg = dxvar_well[nw*2 + w]; +// const double dFo = - dFw - dFg; + //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; +// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit; +// step = std::min(step, 1.0); +// std::cout << step << std::endl; +// F[Water] = xvar_well_old[nw + w] - step*dFw; +// F[Gas] = xvar_well_old[2*nw + w] - step*dFg; +// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo; +// double sumF = F[Water]+F[Gas]+F[Oil]; +// F[Water] /= sumF; +// F[Gas] /= sumF; +// F[Oil] /= sumF; +// well_state.wellSolutions()[nw + w] = F[Water]; +// well_state.wellSolutions()[2 * nw + w] = F[Gas]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + //const int sign1 = dxvar_well[w] > 0 ? 1: -1; + //const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dTotalRateLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dxvar_well[w]; + + switch (wells().type[w]) { + case INJECTOR: + for (int p = 0; p < np; ++p) { + const double comp_frac = wells().comp_frac[np*w + p]; + //if (comp_frac > 0) { + well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w]; + //} + + } + break; + case PRODUCER: + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p]; + } + break; + } + } + break; + case SURFACE_RATE: + { + const int sign1 = dxvar_well[w] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; + //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; + //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; + well_state.bhp()[w] = well_state.wellSolutions()[w]; + + if (wells().type[w]==PRODUCER) { + + double F_target = 0.0; + for (int p = 0; p < np; ++p) { + F_target += wells().comp_frac[np*w + p] * F[p]; + } + for (int p = 0; p < np; ++p) { + //std::cout << F[p] << std::endl; + well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target; + } + } else { + + for (int p = 0; p < np; ++p) { + //std::cout << wells().comp_frac[np*w + p] << " " < 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dBHPLimit); + well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited; + //const int sign = (dxvar_well1[w] < 0) ? -1 : 1; + //well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ; + well_state.bhp()[w] = well_state.wellSolutions()[w]; + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate; + } + } + break; + + } + } + + + + + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + //Loop over all wells +#pragma omp parallel for schedule(static) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + const int nwc = well_controls_get_num(wc); + //Loop over all controls until we find a THP control + //that specifies what we need... + //Will only update THP for wells with THP control + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(wc, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + double alq = well_controls_iget_alq(wc, ctrl_index); + int table_id = well_controls_iget_vfp(wc, ctrl_index); + + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp); + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + //Assume only one THP control specified for each well + break; + } + } + } + } + } + + + + + + template + void + StandardWellsDense:: + updateWellControls(WellState& xw) const + { + if( !localWellsActive() ) return ; + + std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; +#pragma omp parallel for schedule(dynamic) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, wells().type[w], wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + // We disregard terminal_ouput here as with it only messages + // for wells on one process will be printed. + std::ostringstream ss; + ss << "Switching control mode for well " << wells().name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + OpmLog::info(ss.str()); + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + } + + // Updating well state and primary variables. + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + } + } + } else if (well_type == PRODUCER) { + + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + + break; + } + std::vector g = {1,1,0.01}; + if (true) { + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + const WellType& well_type = wells().type[w]; + xw.wellSolutions()[w] = 0.0; + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; + } + } else { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; + } + } + } + break; + + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + xw.wellSolutions()[w] = xw.bhp()[w]; + + } + break; + } + } +// const WellType& well_type = wells().type[w]; +// double tot_well_rate = 0.0; +// for (int p = 0; p < np; ++p) { +// tot_well_rate += g[p] * xw.wellRates()[np*w + p]; +// } +// if (well_type == INJECTOR && tot_well_rate == 0) { +// xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; +// xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; +// } + + } + + } + + + + + + template + void + StandardWellsDense:: + addWellFluxEq(const std::vector& cq_s, + const SolutionState& state, + const double dt, + LinearisedBlackoilResidual& residual) + { + if( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + + double volume = 0.002831684659200; // 0.1 cu ft; + const Vector vol_dt = Vector::Constant(nw,volume/dt); + std::vector F = wellVolumeFractions(state); + //std::cout << F0_[0] << std::endl; + //std::cout << F[0] << std::endl; + + ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); + } + + residual.well_flux_eq = qs; + //std::cout << qs.value() << std::endl; + } + + + + + + template + void + StandardWellsDense::addWellControlEq(const SolutionState& state, + const WellState& xw, + const Vector& aliveWells, + LinearisedBlackoilResidual& residual) + { + if( ! localWellsActive() ) return; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + ADB aqua = ADB::constant(Vector::Zero(nw)); + ADB liquid = ADB::constant(Vector::Zero(nw)); + ADB vapour = ADB::constant(Vector::Zero(nw)); + + if ((*active_)[Water]) { + aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); + } + if ((*active_)[Oil]) { + liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); + } + if ((*active_)[Gas]) { + vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); + } + + //THP calculation variables + std::vector inj_table_id(nw, -1); + std::vector prod_table_id(nw, -1); + Vector thp_inj_target_v = Vector::Zero(nw); + Vector thp_prod_target_v = Vector::Zero(nw); + Vector alq_v = Vector::Zero(nw); + + //Hydrostatic correction variables + Vector rho_v = Vector::Zero(nw); + Vector vfp_ref_depth_v = Vector::Zero(nw); + + //Target vars + Vector bhp_targets = Vector::Zero(nw); + Vector rate_targets = Vector::Zero(nw); + Eigen::SparseMatrix rate_distr(nw, np*nw); + + //Selection variables + std::vector bhp_elems; + std::vector thp_inj_elems; + std::vector thp_prod_elems; + std::vector rate_elems; + + //Run through all wells to calculate BHP/RATE targets + //and gather info about current control + for (int w = 0; w < nw; ++w) { + auto wc = wells().ctrls[w]; + + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + bhp_elems.push_back(w); + bhp_targets(w) = well_controls_iget_target(wc, current); + rate_targets(w) = -1e100; + } + break; + + case THP: + { + const int perf = wells().well_connpos[w]; + rho_v[w] = wellPerforationDensities()[perf]; + + const int table_id = well_controls_iget_vfp(wc, current); + const double target = well_controls_iget_target(wc, current); + + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + inj_table_id[w] = table_id; + thp_inj_target_v[w] = target; + alq_v[w] = -1e100; + + vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + + thp_inj_elems.push_back(w); + } + else if (well_type == PRODUCER) { + prod_table_id[w] = table_id; + thp_prod_target_v[w] = target; + alq_v[w] = well_controls_iget_alq(wc, current); + + vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + + thp_prod_elems.push_back(w); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well"); + } + bhp_targets(w) = -1e100; + rate_targets(w) = -1e100; + } + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + rate_elems.push_back(w); + // RESERVOIR and SURFACE rates look the same, from a + // high-level point of view, in the system of + // simultaneous linear equations. + + const double* const distr = + well_controls_iget_distr(wc, current); + + for (int p = 0; p < np; ++p) { + rate_distr.insert(w, p*nw + w) = distr[p]; + } + + bhp_targets(w) = -1.0e100; + rate_targets(w) = well_controls_iget_target(wc, current); + } + break; + } + } + + //Calculate BHP target from THP + const ADB thp_inj_target = ADB::constant(thp_inj_target_v); + const ADB thp_prod_target = ADB::constant(thp_prod_target_v); + const ADB alq = ADB::constant(alq_v); + const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target); + const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq); + + //Perform hydrostatic correction to computed targets + const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity_); + const ADB dp = ADB::constant(dp_v); + const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); + const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); + + //Calculate residuals + const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj; + const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod; + const ADB bhp_residual = state.bhp - bhp_targets; + const ADB rate_residual = rate_distr * state.qs - rate_targets; + + //Select the right residual for each well + residual.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) + + superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) + + superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) + + superset(subset(rate_residual, rate_elems), rate_elems, nw); + + // For wells that are dead (not flowing), and therefore not communicating + // with the reservoir, we set the equation to be equal to the well's total + // flow. This will be a solution only if the target rate is also zero. + Eigen::SparseMatrix rate_summer(nw, np*nw); + for (int w = 0; w < nw; ++w) { + for (int phase = 0; phase < np; ++phase) { + rate_summer.insert(w, phase*nw + w) = 1.0; + } + } + Selector alive_selector(aliveWells, Selector::NotEqualZero); + residual.well_eq = alive_selector.select(residual.well_eq, rate_summer * state.qs); + // OPM_AD_DUMP(residual_.well_eq); + } + + + + + + template + void + StandardWellsDense::computeWellPotentials(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state0, + WellState& well_state) + { + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + Vector bhps = Vector::Zero(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* ctrl = wells().ctrls[w]; + const int nwc = well_controls_get_num(ctrl); + //Loop over all controls until we find a BHP control + //or a THP control that specifies what we need. + //Pick the value that gives the most restrictive flow + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + + if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { + bhps[w] = well_controls_iget_target(ctrl, ctrl_index); + } + + if (well_controls_iget_type(ctrl, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(ctrl, ctrl_index); + const double& thp = well_controls_iget_target(ctrl, ctrl_index); + const double& alq = well_controls_iget_alq(ctrl, ctrl_index); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + // apply the strictest of the bhp controlls i.e. smallest bhp for injectors + if ( bhp < bhps[w]) { + bhps[w] = bhp; + } + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + // apply the strictest of the bhp controlls i.e. largest bhp for producers + if ( bhp > bhps[w]) { + bhps[w] = bhp; + } + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + } + } + + } + + // use bhp limit from control + state0.bhp = ADB::constant(bhps); + + // compute well potentials + Vector aliveWells; + std::vector well_potentials; + computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_potentials); + + // store well potentials in the well state + // transform to a single vector instead of separate vectors pr phase + const int nperf = wells().well_connpos[nw]; + Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np); + } + well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); + } + + + + + + void + StandardWellsDense::variableStateWellIndices(std::vector& indices, + int& next) const + { + indices[Qs] = next++; + //indices[Bhp] = next++; + } + + + + + template + void + StandardWellsDense:: + variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state, + WellState &xw) const + { + + state.wellVariables = std::move(vars[indices[Qs]]); + //std::cout << "state.wellVariables " << state.wellVariables.value() << std::endl; + + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + Vector isBHPControlled = Vector::Zero(nw); + Vector isReservoirRateControlled = Vector::Zero(nw); + Vector isSurfaceRateControlled = Vector::Zero(nw); + + Vector isInjector = Vector::Zero(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + + switch (well_controls_iget_type(wc, current)) { + + case BHP: + isBHPControlled[w] = 1; + break; + + case RESERVOIR_RATE: + isReservoirRateControlled[w] = 1; + + break; + case SURFACE_RATE: + isSurfaceRateControlled[w] = 1; + break; + } + + if (wells().type[w] == INJECTOR) { + isInjector[w] = 1; + } + } + + Vector ones = Vector::Constant(nw,1); + + + const ADB& xvar_well1 = subset(state.wellVariables,Span(nw,1,0)); + //std::cout << xvar_well1.value() << std::endl; + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + + state.bhp = isBHPControlled * bhp + (1-isBHPControlled) * xvar_well1; + + const ADB Fw = subset(state.wellVariables,Span(nw,1,nw)); + const ADB Fg = subset(state.wellVariables,Span(nw,1,nw*2)); + const ADB Fo = ones - Fw - Fg; + + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + + V target_rates = V::Zero(nw); + V isNoFlow = V::Zero(nw); + std::vector distri(np); + for (int p = 0; p < np; ++p) { + distri[p] = V::Zero(nw); + } + + + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + const double* distr = well_controls_iget_distr(wc, current); + + + target_rates[w] = well_controls_iget_target(wc, current); + if (target_rates[w]==0) + isNoFlow[w] = 1; + for (int p = 0; p < np; ++p) { + const double comp_frac = wells().comp_frac[np*w + p]; + distri[p][w] = distr[p]; + //std::cout << "distr " << distr[p] << " comp_frac " << comp_frac << std::endl; + } + + } + + std::vector F(np, ADB::null()); + F[Water] = Fw; + F[Oil] = Fo; + F[Gas] = Fg; // / Vector::Constant(nw,0.01); + + std::vector g(np,Vector::Constant(nw,1.0)); + g[Gas] = Vector::Constant(nw,0.01); +// g[Water] = Vector::Constant(nw,1.0); +// g[Oil] = Vector::Constant(nw,1.0); + for (int p = 0; p < np; ++p) { + const V tmp = (ones - isReservoirRateControlled) * g[p] + isReservoirRateControlled * distri[p]; + F[p] = F[p] / tmp; + } + + ADB targetF = compi.col(pu.phase_pos[0])*F[0]; + for (int p = 1; p < np; ++p) { + targetF += compi.col(pu.phase_pos[p])*F[p]; + } + + + ADB sumF = distri[0] * F[0]; + for (int p = 1; p < np; ++p) { + sumF += distri[p] * (F[p]); + } + Selector noTargetFSelector(targetF.value(), Selector::Zero); + Selector noSumFSelector(sumF.value(), Selector::Zero); + + + //std::cout << target_rates << std::endl; + ADB qs = ADB::constant(V::Zero(nw*np)); + for (int p = 0; p < np; ++p) { + + const int pos = pu.phase_pos[p]; + ADB Q = isSurfaceRateControlled * ADB::constant(compi.col(pos) * target_rates); + //DUMPVAL(Q); + Q += isBHPControlled * isInjector * xvar_well1 * ADB::constant(compi.col(pos)); + //DUMPVAL(Q); + Q += isBHPControlled * (ones - isInjector) * xvar_well1 * F[p]; + //DUMPVAL(Q); + Q += (ones - isInjector) * isSurfaceRateControlled * noTargetFSelector.select( ADB::constant(V::Zero(nw)), (ones - compi.col(pos)) * target_rates * F[p] / targetF); + + //Q += isNoFlow * -1e-24 * F[p]; + + //Q += isSurfaceRateControlled * isInjector * -1e-12 * F[p] * (ones - compi.col(pos)); + + Q += isReservoirRateControlled * F[p] * target_rates; + //DUMPVAL(Q); + qs += superset(Q, Span(nw, 1, p*nw), nw*np); + } + state.qs = qs; + //if (isNoFlow.sum() > 0) + //std::cout << state.qs << std::endl; + + // Bhp. + //state.bhp = + //state.qs = + } + + + + + + std::vector + StandardWellsDense::variableWellStateIndices() const + { + // Black oil model standard is 5 equation. + // For the pure well solve, only the well equations are picked. + std::vector indices(5, -1); + int next = 0; + + variableStateWellIndices(indices, next); + + assert(next == 1); + return indices; + } + + + + + + template + void + StandardWellsDense::variableWellStateInitials(const WellState& xw, + std::vector& vars0) const + { + // Initial well rates. + if ( localWellsActive() ) + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + + // The transpose() below switches the ordering. + //const DataBlock wrates = Eigen::Map(& xw.wellSolutions()[0], nw, np).transpose(); + //const Vector qs = Eigen::Map(wrates.data(), nw*np); + const Vector qs = Eigen::Map(& xw.wellSolutions()[0], xw.wellSolutions().size()); + vars0.push_back(qs); + + // Initial well bottom-hole pressure. +// assert (not xw.bhp().empty()); +// const Vector bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); +// vars0.push_back(bhp); + + } + else + { + // push null states for qs and bhp + vars0.push_back(V()); + //vars0.push_back(V()); + } + } + + template + std::vector> StandardWellsDense::wellVolumeFractions(const SolutionState& state) const + { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + std::vector F(np,ADB::null()); + F[Water] = subset(state.wellVariables,Span(nw,1,nw)); + F[Gas] = subset(state.wellVariables,Span(nw,1,2 * nw)); + F[Oil] = Vector::Constant(nw,1.0) - F[Water] - F[Gas]; + return F; + } + + + template + void StandardWellsDense::computeAccumWells(const SolutionState& state) { + F0_ = wellVolumeFractions(state); + } + + + + void + StandardWellsDense::setStoreWellPerforationFluxesFlag(const bool store_fluxes) + { + store_well_perforation_fluxes_ = store_fluxes; + } + + + + + + const StandardWellsDense::Vector& + StandardWellsDense::getStoredWellPerforationFluxes() const + { + assert(store_well_perforation_fluxes_); + return well_perforation_fluxes_; + } + + + + + + + template + void + StandardWellsDense:: + updateListEconLimited(ScheduleConstPtr schedule, + const int current_step, + const Wells* wells_struct, + const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const + { + const int nw = wells_struct->number_of_wells; + + for (int w = 0; w < nw; ++w) { + // flag to check if the mim oil/gas rate limit is violated + bool rate_limit_violated = false; + const std::string& well_name = wells_struct->name[w]; + const Well* well_ecl = schedule->getWell(well_name); + const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step); + + // economic limits only apply for production wells. + if (wells_struct->type[w] != PRODUCER) { + continue; + } + + // if no limit is effective here, then continue to the next well + if ( !econ_production_limits.onAnyEffectiveLimit() ) { + continue; + } + // for the moment, we only handle rate limits, not handling potential limits + // the potential limits should not be difficult to add + const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); + if (quantity_limit == WellEcon::POTN) { + const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") + + std::string("All the limits will be evaluated based on RATE. "); + OpmLog::warning("NOT_SUPPORTING_POTN", msg); + } + + const WellMapType& well_map = well_state.wellMap(); + const typename WellMapType::const_iterator i_well = well_map.find(well_name); + assert(i_well != well_map.end()); // should always be found? + const WellMapEntryType& map_entry = i_well->second; + const int well_number = map_entry[0]; + + if (econ_production_limits.onAnyRateLimit()) { + rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number); + } + + if (rate_limit_violated) { + if (econ_production_limits.endRun()) { + const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") + + std::string("the program will keep running after ") + well_name + std::string(" is closed"); + OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); + } + + if (econ_production_limits.validFollowonWell()) { + OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); + } + + if (well_ecl->getAutomaticShutIn()) { + list_econ_limited.addShutWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); + OpmLog::info(msg); + } else { + list_econ_limited.addStoppedWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); + OpmLog::info(msg); + } + // the well is closed, not need to check other limits + continue; + } + + // checking for ratio related limits, mostly all kinds of ratio. + bool ratio_limits_violated = false; + RatioCheckTuple ratio_check_return; + + if (econ_production_limits.onAnyRatioLimit()) { + ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry); + ratio_limits_violated = std::get<0>(ratio_check_return); + } + + if (ratio_limits_violated) { + const bool last_connection = std::get<1>(ratio_check_return); + const int worst_offending_connection = std::get<2>(ratio_check_return); + + const int perf_start = map_entry[1]; + + assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2])); + + const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection]; + list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); + const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") + + well_name + std::string(" will be closed due to economic limit"); + OpmLog::info(msg); + + if (last_connection) { + list_econ_limited.addShutWell(well_name); + const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); + OpmLog::info(msg2); + } + } + + } + } + + + + + + template + bool + StandardWellsDense:: + checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const int well_number) const + { + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int np = well_state.numPhases(); + + if (econ_production_limits.onMinOilRate()) { + assert((*active_)[Oil]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double min_oil_rate = econ_production_limits.minOilRate(); + if (std::abs(oil_rate) < min_oil_rate) { + return true; + } + } + + if (econ_production_limits.onMinGasRate() ) { + assert((*active_)[Gas]); + const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; + const double min_gas_rate = econ_production_limits.minGasRate(); + if (std::abs(gas_rate) < min_gas_rate) { + return true; + } + } + + if (econ_production_limits.onMinLiquidRate() ) { + assert((*active_)[Oil]); + assert((*active_)[Water]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + const double min_liquid_rate = econ_production_limits.minLiquidRate(); + if (std::abs(liquid_rate) < min_liquid_rate) { + return true; + } + } + + if (econ_production_limits.onMinReservoirFluidRate()) { + OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); + } + + return false; + } + + + + + + template + StandardWellsDense::RatioCheckTuple + StandardWellsDense:: + checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const + { + // TODO: not sure how to define the worst-offending connection when more than one + // ratio related limit is violated. + // The defintion used here is that we define the violation extent based on the + // ratio between the value and the corresponding limit. + // For each violated limit, we decide the worst-offending connection separately. + // Among the worst-offending connections, we use the one has the biggest violation + // extent. + + bool any_limit_violated = false; + bool last_connection = false; + int worst_offending_connection = INVALIDCONNECTION; + double violation_extent = -1.0; + + if (econ_production_limits.onMaxWaterCut()) { + const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry); + bool water_cut_violated = std::get<0>(water_cut_return); + if (water_cut_violated) { + any_limit_violated = true; + const double violation_extent_water_cut = std::get<3>(water_cut_return); + if (violation_extent_water_cut > violation_extent) { + violation_extent = violation_extent_water_cut; + worst_offending_connection = std::get<2>(water_cut_return); + last_connection = std::get<1>(water_cut_return); + } + } + } + + if (econ_production_limits.onMaxGasOilRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxWaterGasRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxGasLiquidRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); + } + + if (any_limit_violated) { + assert(worst_offending_connection >=0); + assert(violation_extent > 1.); + } + + return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + + + + + template + StandardWellsDense::RatioCheckTuple + StandardWellsDense:: + checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const + { + bool water_cut_limit_violated = false; + int worst_offending_connection = INVALIDCONNECTION; + bool last_connection = false; + double violation_extent = -1.0; + + const int np = well_state.numPhases(); + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int well_number = map_entry[0]; + + assert((*active_)[Oil]); + assert((*active_)[Water]); + + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + double water_cut; + if (std::abs(liquid_rate) != 0.) { + water_cut = water_rate / liquid_rate; + } else { + water_cut = 0.0; + } + + const double max_water_cut_limit = econ_production_limits.maxWaterCut(); + if (water_cut > max_water_cut_limit) { + water_cut_limit_violated = true; + } + + if (water_cut_limit_violated) { + // need to handle the worst_offending_connection + const int perf_start = map_entry[1]; + const int perf_number = map_entry[2]; + + std::vector water_cut_perf(perf_number); + for (int perf = 0; perf < perf_number; ++perf) { + const int i_perf = perf_start + perf; + const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; + const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; + const double liquid_perf_rate = oil_perf_rate + water_perf_rate; + if (std::abs(liquid_perf_rate) != 0.) { + water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; + } else { + water_cut_perf[perf] = 0.; + } + } + + last_connection = (perf_number == 1); + if (last_connection) { + worst_offending_connection = 0; + violation_extent = water_cut_perf[0] / max_water_cut_limit; + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + double max_water_cut_perf = 0.; + for (int perf = 0; perf < perf_number; ++perf) { + if (water_cut_perf[perf] > max_water_cut_perf) { + worst_offending_connection = perf; + max_water_cut_perf = water_cut_perf[perf]; + } + } + + assert(max_water_cut_perf != 0.); + assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); + + violation_extent = max_water_cut_perf / max_water_cut_limit; + } + + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + +} // namespace Opm diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index ac7fbf428..b796e80d9 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -105,12 +106,67 @@ namespace Opm well_potentials_.clear(); well_potentials_.resize(nperf * np, 0.0); + well_solutions_.clear(); + well_solutions_.resize(nw * np, 0.0); + std::vector g = {1.0,1.0,0.01}; + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells->ctrls[w]; + + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = current_controls_[w]; + const WellType& well_type = wells->type[w]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + well_solutions_[w] += wellRates()[np*w + p] * wells->comp_frac[np*w + p]; + } + } else { + for (int p = 0; p < np; ++p) { + well_solutions_[w] += g[p] * wellRates()[np*w + p]; + } + } + } + break; + + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + wellSolutions()[w] = bhp()[w]; + + } + break; + } + assert(np == 3); + double total_rates = 0.0; + for (int p = 0; p < np; ++p) { + total_rates += g[p] * wellRates()[np*w + p]; + } + + //if(std::abs(total_rates) > 0) { + // wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water; + // wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas + //} else { + wellSolutions()[nw + w] = wells->comp_frac[np*w + Water]; + wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas]; + //} + + + } + + // intialize wells that have been there before // order may change so the mapping is based on the well name if( ! prevState.wellMap().empty() ) { typedef typename WellMapType :: const_iterator const_iterator; const_iterator end = prevState.wellMap().end(); + int nw_old = prevState.bhp().size(); for (int w = 0; w < nw; ++w) { std::string name( wells->name[ w ] ); const_iterator it = prevState.wellMap().find( name ); @@ -123,11 +179,31 @@ namespace Opm bhp()[ newIndex ] = prevState.bhp()[ oldIndex ]; // wellrates + double total_well_rates = 0.0; + for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i 0) { + for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; itype[w] == PRODUCER && std::abs(total_well_rates) > 0.0) { + for( int i=0; i& wellPotentials() { return well_potentials_; } const std::vector& wellPotentials() const { return well_potentials_; } + /// One rate per phase and well connection. + std::vector& wellSolutions() { return well_solutions_; } + const std::vector& wellSolutions() const { return well_solutions_; } + data::Wells report() const override { data::Wells res = WellState::report(); @@ -251,6 +332,7 @@ namespace Opm std::vector perfphaserates_; std::vector current_controls_; std::vector well_potentials_; + std::vector well_solutions_; }; } // namespace Opm From 17541817616add990c6297027ead4ef11e2715d8 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Aug 2016 14:31:25 +0200 Subject: [PATCH 2/8] BUGFIX Make the wellVariables constant The wellVariables where left out in makeConstantState() With this fix norne runs --- opm/autodiff/BlackoilModelEbos.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 6ce28724a..f2a0f0816 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -951,6 +951,7 @@ namespace Opm { ADB& pp = state.canonical_phase_pressures[canphase]; pp = ADB::constant(pp.value()); } + state.wellVariables = ADB::constant(state.wellVariables.value()); } void setupLegacyState(SolutionState& state, From dec60a8bd8c88081b492e3e8214fb98613547cbc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Aug 2016 14:34:27 +0200 Subject: [PATCH 3/8] The well model no uses the fluidState and fluidsystem from ebos The well model is modified to use the fluidState and fluidsystem from ebos. UpdateLegacyState() is no longer needed. --- opm/autodiff/BlackoilModelEbos.hpp | 409 +++++++++++++++++++++-- opm/autodiff/StandardWellsDense.hpp | 27 +- opm/autodiff/StandardWellsDense_impl.hpp | 60 ++-- 3 files changed, 433 insertions(+), 63 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index f2a0f0816..01aa91268 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -300,7 +300,7 @@ namespace Opm { makeConstantState(state0); // Compute initial accumulation contributions // and well connection pressures. - wellModel().computeWellConnectionPressures(state0, well_state); + computeWellConnectionPressures(state0, well_state); wellModel().computeAccumWells(state0); } @@ -309,17 +309,18 @@ namespace Opm { return iter_report; } + double dt = timer.currentStepLength(); std::vector mob_perfcells; std::vector b_perfcells; - wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); + //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && iterationIdx == 0 ) { // solve the well equations as a pre-processing step iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, state, well_state); } V aliveWells; - std::vector cq_s; - wellModel().computeWellFluxDense(state, mob_perfcells, b_perfcells, cq_s); + std::vector cq_s; + computeWellFluxDense(state, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); wellModel().addWellFluxEq(cq_s, state, dt, residual_); addWellContributionToMassBalanceEq(cq_s, state, well_state); @@ -333,6 +334,341 @@ namespace Opm { return iter_report; } + typedef DenseAd::Evaluation Eval; + EvalWell extendEval(Eval in) const { + EvalWell out = 0.0; + out.value = in.value; + for(int i = 0;i<3;++i) { + out.derivatives[i] = in.derivatives[flowToEbosPvIdx(i)]; + } + return out; + } + + template + void + computeWellFluxDense(const SolutionState& state, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + V Tw = Eigen::Map(wellModel().wells().WI, nperf); + const std::vector& well_cells = wellModel().wellOps().well_cells; + std::vector well_id(nperf); + + // pressure diffs computed already (once per step, not changing per iteration) + const V& cdp = wellModel().wellPerforationPressureDiffs(); + + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + std::vector cmix_s_ADB = wellModel().wellVolumeFractions(state); + + for (int w = 0; w < nw; ++w) { + + EvalWell bhp = wellModel().extractDenseADWell(state.bhp,w); + + + // TODO: fix for 2-phase case + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + cmix_s[phase] = wellModel().extractDenseADWell(cmix_s_ADB[phase],w); + } + + //std::cout <<"cmix gas "<< w<< " "< b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + } + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + cdp[perf]; + EvalWell drawdown = pressure - well_pressure; + + // injection perforations + if ( drawdown.value > 0 ) { + + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) + continue; + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + } + + + if ((active_)[Oil] && (active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs * cq_psOil; + cq_ps[oilpos] += rv * cq_psGas; + } + + + // map to ADB + for (int phase = 0; phase < np; ++phase) { + cq_s_dense[phase][perf] = cq_ps[phase]; + } + + } else { + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) + continue; + + // Using total mobilities + EvalWell total_mob_dense = mob_perfcells_dense[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob_dense += mob_perfcells_dense[phase]; + } + // injection perforations total volume rates + const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio = 0.0; + if ((active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if ((active_)[Oil] && (active_)[Gas]) { + EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); + EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + EvalWell rvPerf = 0.0; + if (cmix_s[gaspos] > 0) + rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; + + if (rvPerf.value > rvSatEval.value) { + rvPerf = rvSatEval; + //rvPerf.value = rvSatEval.value; + } + + EvalWell rsPerf = 0.0; + if (cmix_s[oilpos] > 0) + rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; + + if (rsPerf.value > rsSatEval.value) { + //rsPerf = 0.0; + rsPerf= rsSatEval; + } + + // Incorporate RS/RV factors if both oil and gas active + const EvalWell d = 1.0 - rvPerf * rsPerf; + + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < cq_s2; + //Vector aliveWells; + //computeWellFlux(state,mob_perfcells,b_perfcells, aliveWells,cq_s2); + + //for (int phase = 0; phase < np; ++phase) { + //if( !(((cq_s[phase].value() - cq_s2[phase].value()).abs()<1e-10).all()) ) { + //std::cout << "phase " << phase << std::endl; + //std::cout << cq_s2[phase].value() << std::endl; + //std::cout << cq_s[phase].value() << std::endl; + //} + //} + } + void + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + 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, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); + const V& pdepth = subset(depth, wellModel().wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + wellModel().computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + template + void + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const std::vector& well_cells = wellModel().wellOps().well_cells; + const PhaseUsage& pu = fluid_.phaseUsage(); + const int np = fluid_.numPhases(); + b_perf.resize(nperf*np); + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = phaseCondition_[well_cells[perf]]; + } + + // Compute the average pressure in each well block + const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); + //V avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = well_cells[perf]; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeOil()) { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rv = fs.Rv().value; + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeGas()) { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rs = fs.Rs().value; + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } + } + + +// // Use cell values for the temperature as the wells don't knows its temperature yet. +// const ADB perf_temp = subset(state.temperature, well_cells); + +// // Compute b, rsmax, rvmax values for perforations. +// // Evaluate the properties using average well block pressures +// // and cell values for rs, rv, phase condition and temperature. +// const ADB avg_press_ad = ADB::constant(avg_press); +// // const std::vector& pc = phaseCondition(); + +// DataBlock b(nperf, pu.num_phases); +// if (pu.phase_used[BlackoilPhases::Aqua]) { +// const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; +// } +// assert((*active_)[Oil]); +// const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); +// if (pu.phase_used[BlackoilPhases::Liquid]) { +// const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null(); +// const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; +// } +// if (pu.phase_used[BlackoilPhases::Vapour]) { +// const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null(); +// const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; +// } +// if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { +// const V rssat = fluid_.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); +// //rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + +// const V rvsat = fluid_.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); +// //rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); +// } + +// // b is row major, so can just copy data. +// //b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + + /// \brief Compute the residual norms of the mass balance for each phase, /// the well flux, and the well equation. @@ -710,8 +1046,23 @@ namespace Opm { for ( int idx = 0; idx < np; ++idx ) { - const ADB& tempB = rq_[idx].b; - B.col(idx) = 1./tempB.value(); + V b(nc); + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + } + B.col(idx) = b; + } + + + for ( int idx = 0; idx < np; ++idx ) + { + //const ADB& tempB = rq_[idx].b; + //B.col(idx) = 1./tempB.value(); R.col(idx) = residual_.material_balance_eq[idx].value(); tempV.col(idx) = R.col(idx).abs()/pv; } @@ -1460,7 +1811,7 @@ namespace Opm { std::move(pAdbJacs[phaseIdx])); rq_[phaseIdx].mob = ADB::function(std::move(mobVal[phaseIdx]), - std::move(mobAdbJacs[phaseIdx])); + std::move(mobAdbJacs[phaseIdx])); rq_[phaseIdx].b = ADB::function(std::move(bVal[phaseIdx]), std::move(bAdbJacs[phaseIdx])); @@ -1522,7 +1873,7 @@ namespace Opm { prevEpisodeIdx = ebosSimulator_.episodeIndex(); convertResults(ebosSimulator_, /*sparsityPattern=*/state.saturation[0]); - updateLegacyState(ebosSimulator_, state); + //updateLegacyState(ebosSimulator_, state); if (param_.update_equations_scaling_) { updateEquationsScaling(); @@ -1536,7 +1887,6 @@ namespace Opm { SolutionState& state, WellState& well_state) { - V aliveWells; const int np = wells().number_of_phases; std::vector cq_s(np, ADB::null()); std::vector indices = wellModel().variableWellStateIndices(); @@ -1550,10 +1900,10 @@ namespace Opm { if (localWellsActive() ){ // If there are non well in the sudomain of the process // thene mob_perfcells_const and b_perfcells_const would be empty - for (int phase = 0; phase < np; ++phase) { - mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); - b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); - } +// for (int phase = 0; phase < np; ++phase) { +// mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); +// b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); +// } } int it = 0; @@ -1567,7 +1917,7 @@ namespace Opm { SolutionState wellSolutionState = state0; wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState, well_state); - wellModel().computeWellFluxDense(wellSolutionState, mob_perfcells_const, b_perfcells_const, cq_s); + computeWellFluxDense(wellSolutionState, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); wellModel().addWellFluxEq(cq_s, wellSolutionState, dt, residual_); //wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); @@ -1672,10 +2022,25 @@ namespace Opm { Eigen::Array B(nc, np); Eigen::Array R(nc, np); Eigen::Array tempV(nc, np); + for ( int idx = 0; idx < np; ++idx ) { - const ADB& tempB = rq_[idx].b; - B.col(idx) = 1./tempB.value(); + V b(nc); + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + } + B.col(idx) = b; + } + + for ( int idx = 0; idx < np; ++idx ) + { + //const ADB& tempB = rq_[idx].b; + //B.col(idx) = 1./tempB.value(); R.col(idx) = residual_.material_balance_eq[idx].value(); tempV.col(idx) = R.col(idx).abs()/pv; } @@ -1835,12 +2200,12 @@ namespace Opm { // TODO: added since the interfaces of the function are different // TODO: for StandardWells and MultisegmentWells - void - computeWellConnectionPressures(const SolutionState& state, - const WellState& well_state) - { - wellModel().computeWellConnectionPressures(state, well_state); - } +// void +// computeWellConnectionPressures(const SolutionState& state, +// const WellState& well_state) +// { +// wellModel().computeWellConnectionPressures(state, well_state); +// } /// \brief Compute the reduction within the convergence check. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 601760fdb..6ffb2c8c8 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -60,11 +60,12 @@ namespace Opm { // --------- Types --------- using ADB = AutoDiffBlock; - typedef DenseAd::Evaluation Eval; + typedef DenseAd::Evaluation EvalWell; //typedef AutoDiffBlock ADB; using Vector = ADB::V; using V = ADB::V; + // copied from BlackoilModelBase // should put to somewhere better using DataBlock = Eigen::Array& b_perfcells, std::vector& cq_s) const; - Eval extractDenseAD(const ADB& data, int i, int j) const; - Eval extractDenseADWell(const ADB& data, int i) const; + EvalWell extractDenseAD(const ADB& data, int i, int j) const; + EvalWell extractDenseADWell(const ADB& data, int i) const; - const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; + const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; template @@ -215,6 +216,15 @@ namespace Opm { template void computeAccumWells(const SolutionState& state); + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav); + protected: bool wells_active_; const Wells* wells_; @@ -246,14 +256,7 @@ namespace Opm { std::vector& rvmax_perf, std::vector& surf_dens_perf); - template - void computeWellConnectionDensitesPressures(const WellState& xw, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf, - const std::vector& depth_perf, - const double grav); + template diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 46fd6b85d..3374fe846 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -397,6 +397,7 @@ namespace Opm std::vector& cq_s) const { if( ! localWellsActive() ) return ; + const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -408,7 +409,7 @@ namespace Opm // pressure diffs computed already (once per step, not changing per iteration) const Vector& cdp = wellPerforationPressureDiffs(); - std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); std::vector cmix_s_ADB = wellVolumeFractions(state); const int oilpos = pu.phase_pos[Oil]; @@ -417,11 +418,12 @@ namespace Opm const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s_ADB[oilpos], well_cells); for (int w = 0; w < nw; ++w) { - Eval bhp = extractDenseADWell(state.bhp,w); + + EvalWell bhp = extractDenseADWell(state.bhp,w); // TODO: fix for 2-phase case - std::vector cmix_s(np,0.0); + std::vector cmix_s(np,0.0); for (int phase = 0; phase < np; ++phase) { cmix_s[phase] = extractDenseADWell(cmix_s_ADB[phase],w); } @@ -431,11 +433,11 @@ namespace Opm for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { const int cell_idx = well_cells[perf]; well_id[perf] = w; - Eval pressure = extractDenseAD(state.pressure, cell_idx, cell_idx); - Eval rs = extractDenseAD(state.rs, cell_idx, cell_idx); - Eval rv = extractDenseAD(state.rv, cell_idx, cell_idx); - std::vector b_perfcells_dense(np, 0.0); - std::vector mob_perfcells_dense(np, 0.0); + EvalWell pressure = extractDenseAD(state.pressure, cell_idx, cell_idx); + EvalWell rs = extractDenseAD(state.rs, cell_idx, cell_idx); + EvalWell rv = extractDenseAD(state.rv, cell_idx, cell_idx); + std::vector b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); for (int phase = 0; phase < np; ++phase) { b_perfcells_dense[phase] = extractDenseAD(b_perfcells[phase], perf, cell_idx); mob_perfcells_dense[phase] = extractDenseAD(mob_perfcells[phase], perf, cell_idx); @@ -443,7 +445,7 @@ namespace Opm } // Pressure drawdown (also used to determine direction of flow) - Eval drawdown = pressure - bhp - cdp[perf]; + EvalWell drawdown = pressure - bhp - cdp[perf]; // injection perforations if ( drawdown.value > 0 ) { @@ -452,9 +454,9 @@ namespace Opm if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) continue; // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, 0.0); + std::vector cq_ps(np, 0.0); for (int phase = 0; phase < np; ++phase) { - const Eval cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); cq_ps[phase] = b_perfcells_dense[phase] * cq_p; } @@ -462,8 +464,8 @@ namespace Opm if ((*active_)[Oil] && (*active_)[Gas]) { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; - const Eval cq_psOil = cq_ps[oilpos]; - const Eval cq_psGas = cq_ps[gaspos]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; cq_ps[gaspos] += rs * cq_psOil; cq_ps[oilpos] += rv * cq_psGas; } @@ -480,15 +482,15 @@ namespace Opm continue; // Using total mobilities - Eval total_mob_dense = mob_perfcells_dense[0]; + EvalWell total_mob_dense = mob_perfcells_dense[0]; for (int phase = 1; phase < np; ++phase) { total_mob_dense += mob_perfcells_dense[phase]; } // injection perforations total volume rates - const Eval cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); // compute volume ratio between connection at standard conditions - Eval volumeRatio = 0.0; + EvalWell volumeRatio = 0.0; if ((*active_)[Water]) { const int watpos = pu.phase_pos[Water]; volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; @@ -498,7 +500,7 @@ namespace Opm const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; - Eval rvPerf = 0.0; + EvalWell rvPerf = 0.0; if (cmix_s[gaspos] > 0) rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; @@ -507,7 +509,7 @@ namespace Opm rvPerf.value = rvSat.value()[w]; } - Eval rsPerf = 0.0; + EvalWell rsPerf = 0.0; if (cmix_s[oilpos] > 0) rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; @@ -517,13 +519,13 @@ namespace Opm } // Incorporate RS/RV factors if both oil and gas active - const Eval d = 1.0 - rvPerf * rsPerf; + const EvalWell d = 1.0 - rvPerf * rsPerf; - const Eval tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; //std::cout << "tmp_oil " < Eval; - Eval + typedef DenseAd::Evaluation EvalWell; + EvalWell StandardWellsDense:: extractDenseAD(const ADB& data, int i, int j) const { - Eval output = 0.0; + EvalWell output = 0.0; output.value = data.value()[i]; const int np = wells().number_of_phases; const std::vector& jac = data.derivative(); @@ -831,12 +833,12 @@ namespace Opm return output; } - typedef DenseAd::Evaluation Eval; - Eval + typedef DenseAd::Evaluation EvalWell; + EvalWell StandardWellsDense:: extractDenseADWell(const ADB& data, int i) const { - Eval output = 0.0; + EvalWell output = 0.0; output.value = data.value()[i]; const int nw = wells().number_of_wells; const int np = wells().number_of_phases; @@ -849,7 +851,7 @@ namespace Opm return output; } - const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const + const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const { typedef typename ADB::M M; const int nLocal = local.size(); From 49f478480d51e6f9089ee492cbdb8b00446f6182 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 17 Aug 2016 09:54:32 +0200 Subject: [PATCH 4/8] Remove usage of state in the model 1) The wellsolution is stored in wellVariables as a vector of DenseAD objects. The wellVariables is computed from the wellstate in setWellVariables() 2) BHP and well fluxes are not stored directly but calculated based on the wellVariables 3) The initial well accumulation term is stored as a vector of doubles. 4) addWellFluxEq() uses DenseAd flux and accumulation terms 5) computePropertiesForWellConnectionPressures() no longer uses state as input. The wellstate is used to get the bhp pressure. 6) The current wellcontrol is set when updated in updateWellControls() in order to use well_controls_get_current_type(wc) insted of well_controls_iget_type(wc, currentIdx) Tested on SPE1, SPE9 and Norne. Effects the convergence but the results are identical. --- opm/autodiff/BlackoilModelEbos.hpp | 428 +++++++++++++++++------ opm/autodiff/StandardWellsDense.hpp | 2 +- opm/autodiff/StandardWellsDense_impl.hpp | 11 +- 3 files changed, 338 insertions(+), 103 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 01aa91268..e80d3dddb 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -172,6 +172,8 @@ namespace Opm { , terminal_output_ (terminal_output) , current_relaxation_(1.0) , isBeginReportStep_(false) + , wellVariables_(wells().number_of_wells * wells().number_of_phases) + , F0_(wells().number_of_wells * wells().number_of_phases) { const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); @@ -286,22 +288,30 @@ namespace Opm { wellModel().updateWellControls(well_state); // Create the primary variables. - SolutionState state(/*numPhases=*/3); - setupLegacyState(state, reservoir_state, well_state); + setWellVariables(well_state); + + //SolutionState state(/*numPhases=*/3); + //setupLegacyState(state, reservoir_state, well_state); // -------- Mass balance equations -------- - assembleMassBalanceEq(timer, iterationIdx, reservoir_state, state); + assembleMassBalanceEq(timer, iterationIdx, reservoir_state); // -------- Well equations ---------- if (iterationIdx == 0) { // Create the (constant, derivativeless) initial state. - SolutionState state0 = state; - makeConstantState(state0); // Compute initial accumulation contributions // and well connection pressures. - computeWellConnectionPressures(state0, well_state); - wellModel().computeAccumWells(state0); + computeWellConnectionPressures(well_state); + + computeAccumWells(); + // Create the (constant, derivativeless) initial state. + //SolutionState state0 = state; + //makeConstantState(state0); + // Compute initial accumulation contributions + // and well connection pressures. + //computeWellConnectionPressures(state0, well_state); + //wellModel().computeAccumWells(state0); } IterationReport iter_report = {false, false, 0, 0}; @@ -314,22 +324,33 @@ namespace Opm { std::vector mob_perfcells; std::vector b_perfcells; //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); - if (param_.solve_welleq_initially_ && iterationIdx == 0 ) { + if (param_.solve_welleq_initially_ && iterationIdx == 0) { // solve the well equations as a pre-processing step - iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, state, well_state); + iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, well_state); } - V aliveWells; std::vector cq_s; - computeWellFluxDense(state, cq_s); - wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); - wellModel().addWellFluxEq(cq_s, state, dt, residual_); - addWellContributionToMassBalanceEq(cq_s, state, well_state); + computeWellFluxDense(cq_s, 4); + updatePerfPhaseRatesAndPressures(cq_s, well_state); + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int p = 0; p < np; ++p) { + for (int w = 0; w < nw; ++w) { + //std::cout << p << " " << w << " " << getQs(w,p).value< + void print(EvalWell in) const { + std::cout << in.value << std::endl; + for (int i = 0; i < in.derivatives.size(); ++i) { + std::cout << in.derivatives[i] << std::endl; + } + } + void - computeWellFluxDense(const SolutionState& state, - std::vector& cq_s) const + setWellVariables(const WellState& xw) { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + wellVariables_[w + nw*phaseIdx] = 0.0; + wellVariables_[w + nw*phaseIdx].value = xw.wellSolutions()[w + nw* phaseIdx]; + wellVariables_[w + nw*phaseIdx].derivatives[np + phaseIdx] = 1.0; + } + } + } + + void + computeAccumWells() { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + F0_[w + nw * phaseIdx] = wellVolumeFraction(w,phaseIdx).value; + } + } + } + + + template + void + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + WellState& xw) const + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + + V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); + + // Update the perforation pressures. + const V& cdp = wellModel().wellPerforationPressureDiffs(); + for (int w = 0; w < nw; ++w ) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + xw.perfPress()[perf] = cdp[perf] + xw.bhp()[w]; + } + } + + } + + void + addWellFluxEq(std::vector cq_s, + const double dt, + const int numBlocks) + { + if( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + double volume = 0.002831684659200; // 0.1 cu ft; + //std::vector F = wellVolumeFractions(state); + //std::cout << F0_[0] << std::endl; + //std::cout << F[0] << std::endl; + //std::cout << "før Ebos" < res_vec(nw); + for (int w = 0; w < nw; ++w) { + + EvalWell res = (wellVolumeFraction(w, p) - F0_[w + nw*p]) * volume / dt; + res += getQs(w, p); + //for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // res -= cq_s[perf*np + p]; + //} + res_vec[w] = res; + + } + ADB tmp = convertToADBWell(res_vec, numBlocks); + qs += superset(tmp,Span(nw,1,p*nw), nw*np); + } + //std::cout << residual_.well_flux_eq << std::endl; + + + //wellModel().convertToADB(res_vec, well_cells, nc, well_id, nw, numBlocks); + //ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wellModel().wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + //qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); + } + + residual_.well_flux_eq = qs; + + //std::cout << "etter Ebos" << residual_.well_flux_eq << std::endl; + + } + const AutoDiffBlock convertToADBWell(const std::vector& local, const int numVars) const + { + typedef typename ADB::M M; + const int nLocal = local.size(); + typename ADB::V value( nLocal ); + //const int numVars = 5; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + Eigen::SparseMatrix matFlux(nLocal,np*nw); + + for( int i=0; i jacs( numVars ); + if (numVars == 4) { + for( int d=0; d(deri[d]); + //jacs[ d ] = M(mat[d]); + } + + jacs[3] = M(matFlux); + //jacs[4] = M(matBHP); + } + else if (numVars == 1) { + jacs[0] = M(matFlux); + //jacs[1] = M(matBHP); + } + //std::cout << numVars << std::endl; + + return ADB::function( std::move( value ), std::move( jacs )); + } + + + void + computeWellFluxDense(std::vector& cq_s, const int numBlocks) const { if( ! localWellsActive() ) return ; - const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -358,22 +532,30 @@ namespace Opm { V Tw = Eigen::Map(wellModel().wells().WI, nperf); const std::vector& well_cells = wellModel().wellOps().well_cells; std::vector well_id(nperf); + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = wellModel().wellPerforationPressureDiffs(); - std::vector> cq_s_dense(np, std::vector(nperf,0.0)); - std::vector cmix_s_ADB = wellModel().wellVolumeFractions(state); + //std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + //std::vector cmix_s_ADB = wellModel().wellVolumeFractions(state); for (int w = 0; w < nw; ++w) { - EvalWell bhp = wellModel().extractDenseADWell(state.bhp,w); + EvalWell bhp = getBhp(w); //wellModel().extractDenseADWell(state.bhp,w); +// std::cout << "well " << w << std::endl; +// std::cout << "bhpF " << std::endl; +// print(bhp); +// std::cout << "state.bhp " << std::endl; +// print(wellModel().extractDenseADWell(state.bhp,w)); // TODO: fix for 2-phase case std::vector cmix_s(np,0.0); for (int phase = 0; phase < np; ++phase) { - cmix_s[phase] = wellModel().extractDenseADWell(cmix_s_ADB[phase],w); + //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + cmix_s[phase] = wellVolumeFraction(w,phase); } //std::cout <<"cmix gas "<< w<< " "< rsmax_perf; std::vector rvmax_perf; std::vector surf_dens_perf; - computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + computePropertiesForWellConnectionPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); @@ -546,10 +728,9 @@ namespace Opm { } - template + template void - computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, + computePropertiesForWellConnectionPressures(const WellState& xw, std::vector& b_perf, std::vector& rsmax_perf, std::vector& rvmax_perf, @@ -578,7 +759,7 @@ namespace Opm { const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_above = perf == wells().well_connpos[w] ? xw.bhp()[w] : perf_press[perf - 1]; const double p_avg = (perf_press[perf] + p_above)/2; double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value; @@ -1264,9 +1445,110 @@ namespace Opm { double current_relaxation_; V dx_old_; + std::vector wellVariables_; + std::vector F0_; + // --------- Protected methods --------- public: + + EvalWell getBhp(const int wellIdx) const { + const WellControls* wc = wells().ctrls[wellIdx]; + if (well_controls_get_current_type(wc) == BHP) { + EvalWell bhp = 0.0; + const double target_rate = well_controls_get_current_target(wc); + bhp.value = target_rate; + return bhp; + } + return wellVariables_[wellIdx]; + } + + EvalWell getQs(const int wellIdx, const int phaseIdx) const { + EvalWell qs = 0.0; + const WellControls* wc = wells().ctrls[wellIdx]; + const int np = fluid_.numPhases(); + const double target_rate = well_controls_get_current_target(wc); +// std::cout << "well info " << std::endl; +// std::cout << wellIdx << " " << wells().type[wellIdx] << " " < g = {1,1,0.01}; + return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); + } + + + + /// return the StandardWells object StandardWellsDense& wellModel() { return well_model_; } const StandardWellsDense& wellModel() const { return well_model_; } @@ -1490,7 +1772,7 @@ namespace Opm { } private: - void convertResults(const Simulator& simulator, const ADB& sparsityPattern) + void convertResults(const Simulator& simulator) { const auto& ebosJac = simulator.model().linearizer().matrix(); const auto& ebosResid = simulator.model().linearizer().residual(); @@ -1579,9 +1861,10 @@ namespace Opm { adbJacs[ flowPhaseIdx ][ pvIdx ].assign( std::move(jacs[ flowPhaseIdx ][ pvIdx ]) ); } // add two "dummy" matrices for the well primary variables + const int numWells = wells().number_of_wells; for( int pvIdx = numPhases; pvIdx < numPhases + 1; ++pvIdx ) { - adbJacs[ flowPhaseIdx ][ pvIdx ] = - sparsityPattern.derivative()[pvIdx]; + adbJacs[ flowPhaseIdx ][ pvIdx ] = ADBJacobianMatrix( numPhases*numWells, numPhases*numWells ); + //sparsityPattern.derivative()[pvIdx]; } } @@ -1832,8 +2115,7 @@ namespace Opm { private: void assembleMassBalanceEq(const SimulatorTimerInterface& timer, const int iterationIdx, - const ReservoirState& reservoirState, - SolutionState& state) + const ReservoirState& reservoirState) { convertInput( iterationIdx, reservoirState, ebosSimulator_ ); @@ -1872,7 +2154,7 @@ namespace Opm { prevEpisodeIdx = ebosSimulator_.episodeIndex(); - convertResults(ebosSimulator_, /*sparsityPattern=*/state.saturation[0]); + convertResults(ebosSimulator_); //updateLegacyState(ebosSimulator_, state); if (param_.update_equations_scaling_) { @@ -1884,43 +2166,22 @@ namespace Opm { IterationReport solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, const double dt, - SolutionState& state, WellState& well_state) { const int np = wells().number_of_phases; + //std::vector cq_s(np, ADB::null()); + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; std::vector cq_s(np, ADB::null()); - std::vector indices = wellModel().variableWellStateIndices(); - SolutionState state0 = state; WellState well_state0 = well_state; - makeConstantState(state0); - - std::vector mob_perfcells_const(np, ADB::null()); - std::vector b_perfcells_const(np, ADB::null()); - - if (localWellsActive() ){ - // If there are non well in the sudomain of the process - // thene mob_perfcells_const and b_perfcells_const would be empty -// for (int phase = 0; phase < np; ++phase) { -// mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); -// b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); -// } - } int it = 0; bool converged; do { // bhp and Q for the wells - std::vector vars0; - vars0.reserve(1); - wellModel().variableWellStateInitials(well_state, vars0); - std::vector vars = ADB::variables(vars0); - - SolutionState wellSolutionState = state0; - wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState, well_state); - computeWellFluxDense(wellSolutionState, cq_s); - wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); - wellModel().addWellFluxEq(cq_s, wellSolutionState, dt, residual_); - //wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); + computeWellFluxDense(cq_s, 1); + updatePerfPhaseRatesAndPressures(cq_s, well_state); + addWellFluxEq(cq_s, dt, 1); converged = getWellConvergence(it); if (converged) { @@ -1945,37 +2206,10 @@ namespace Opm { assert(dx.size() == total_residual_v.size()); wellModel().updateWellState(dx.array(), dpMaxRel(), well_state); wellModel().updateWellControls(well_state); + setWellVariables(well_state); } } while (it < 15); - if (converged) { - OpmLog::note("well converged iter: " + std::to_string(it)); - const int nw = wells().number_of_wells; - { - // We will set the bhp primary variable to the new ones, - // but we do not change the derivatives here. - ADB::V new_bhp = Eigen::Map(well_state.bhp().data(), nw); - // Avoiding the copy below would require a value setter method - // in AutoDiffBlock. - std::vector old_derivs = state.bhp.derivative(); - state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs)); - } - { - // Need to reshuffle well rates, from phase running fastest - // to wells running fastest. - // The transpose() below switches the ordering. - const DataBlock wrates = Eigen::Map(well_state.wellRates().data(), nw, np).transpose(); - ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); - std::vector old_derivs = state.qs.derivative(); - state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); - } - - const ADB::V new_well_var = Eigen::Map(& well_state.wellSolutions()[0], well_state.wellSolutions().size()); - std::vector old_derivs = state.wellVariables.derivative(); - state.wellVariables = ADB::function(std::move(new_well_var), std::move(old_derivs)); - //computeWellConnectionPressures(state, well_state); - } - if (!converged) { well_state = well_state0; } @@ -1987,9 +2221,7 @@ namespace Opm { void - addWellContributionToMassBalanceEq(const std::vector& cq_s, - const SolutionState& state, - const WellState& xw) + addWellContributionToMassBalanceEq(const std::vector& cq_s) { if ( !localWellsActive() ) { @@ -2003,7 +2235,7 @@ namespace Opm { const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { residual_.material_balance_eq[phase] -= superset(cq_s[phase], wellModel().wellOps().well_cells, nc); - } + } } diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 6ffb2c8c8..86ec1034b 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -143,7 +143,7 @@ namespace Opm { WellState& well_state); template - void updateWellControls(WellState& xw) const; + void updateWellControls(WellState& xw); // TODO: should LinearisedBlackoilResidual also be a template class? template diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 3374fe846..a6d4a9646 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -901,7 +901,6 @@ namespace Opm - template void StandardWellsDense:: @@ -920,6 +919,7 @@ namespace Opm const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; + Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); for (int phase = 1; phase < np; ++phase) { cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); @@ -1194,7 +1194,7 @@ namespace Opm template void StandardWellsDense:: - updateWellControls(WellState& xw) const + updateWellControls(WellState& xw) { if( !localWellsActive() ) return ; @@ -1205,7 +1205,7 @@ namespace Opm const int nw = wells().number_of_wells; #pragma omp parallel for schedule(dynamic) for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; + WellControls* wc = wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1242,6 +1242,8 @@ namespace Opm current = xw.currentControls()[w]; } + well_controls_set_current( wc, current); + // Updating well state and primary variables. // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE const double target = well_controls_iget_target(wc, current); @@ -1414,8 +1416,9 @@ namespace Opm qs += superset((F[phase]-F0_[phase]) * vol_dt, Span(nw,1,phase*nw), nw*np); } + residual.well_flux_eq = qs; - //std::cout << qs.value() << std::endl; + //std::cout << "etter dense " < Date: Thu, 18 Aug 2016 08:58:07 +0200 Subject: [PATCH 5/8] Minor convergence improvments - set current control when initializing the wellstate - re calculate wellVariable after well control has changed. --- opm/autodiff/BlackoilModelEbos.hpp | 2 +- opm/autodiff/StandardWellsDense_impl.hpp | 251 +++++++++--------- .../WellStateFullyImplicitBlackoil.hpp | 19 +- 3 files changed, 140 insertions(+), 132 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index e80d3dddb..f72d42bc1 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -287,7 +287,7 @@ namespace Opm { // get reasonable initial conditions for the wells wellModel().updateWellControls(well_state); - // Create the primary variables. + // Set the primary variables for the wells setWellVariables(well_state); //SolutionState state(/*numPhases=*/3); diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index a6d4a9646..cedcf233b 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -74,7 +74,6 @@ namespace Opm - StandardWellsDense::StandardWellsDense(const Wells* wells_arg) : wells_active_(wells_arg!=nullptr) , wells_(wells_arg) @@ -1235,150 +1234,156 @@ namespace Opm // for wells on one process will be printed. std::ostringstream ss; ss << "Switching control mode for well " << wells().name[w] - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; OpmLog::info(ss.str()); xw.currentControls()[w] = ctrl_index; current = xw.currentControls()[w]; - } + well_controls_set_current( wc, current); - well_controls_set_current( wc, current); - // Updating well state and primary variables. - // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE - const double target = well_controls_iget_target(wc, current); - const double* distr = well_controls_iget_distr(wc, current); - switch (well_controls_iget_type(wc, current)) { - case BHP: - xw.bhp()[w] = target; - break; - case THP: { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; + // Updating well state and primary variables if constraint is broken - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; - if ((*active_)[ Water ]) { - aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if ((*active_)[ Oil ]) { - liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if ((*active_)[ Gas ]) { - vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; } - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; + case SURFACE_RATE: + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + //if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + //} + } + } else if (well_type == PRODUCER) { - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } - xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + + break; } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); - - xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - break; - } - - case RESERVOIR_RATE: - // No direct change to any observable quantity at - // surface condition. In this case, use existing - // flow rates as initial conditions as reservoir - // rate acts only in aggregate. - break; - - case SURFACE_RATE: - // assign target value as initial guess for injectors and - // single phase producers (orat, grat, wrat) - const WellType& well_type = wells().type[w]; - if (well_type == INJECTOR) { + std::vector g = {1,1,0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + const double* distr = well_controls_iget_distr(wc, current); for (int phase = 0; phase < np; ++phase) { - const double& compi = wells().comp_frac[np * w + phase]; - if (compi > 0.0) { - xw.wellRates()[np*w + phase] = target * compi; + g[phase] = distr[phase]; + } + } + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + const WellType& well_type = wells().type[w]; + xw.wellSolutions()[w] = 0.0; + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; + } + } else { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; } } - } else if (well_type == PRODUCER) { + } + break; - // only set target as initial rates for single phase - // producers. (orat, grat and wrat, and not lrat) - // lrat will result in numPhasesWithTargetsUnderThisControl == 2 - int numPhasesWithTargetsUnderThisControl = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - numPhasesWithTargetsUnderThisControl += 1; - } - } - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { - xw.wellRates()[np*w + phase] = target * distr[phase]; - } - } + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + xw.wellSolutions()[w] = xw.bhp()[w]; + } + break; + } + + double tot_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + tot_well_rate += g[p] * xw.wellRates()[np*w + p]; + } + if(std::abs(tot_well_rate) > 0) { + xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; //wells->comp_frac[np*w + Water]; // Water; + xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; //wells->comp_frac[np*w + Gas]; //Gas } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - - break; - } - std::vector g = {1,1,0.01}; - if (true) { - switch (well_controls_iget_type(wc, current)) { - case BHP: - { - const WellType& well_type = wells().type[w]; - xw.wellSolutions()[w] = 0.0; - if (well_type == INJECTOR) { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; - } - } else { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; - } + //xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; + //xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; } } - break; - - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - { - xw.wellSolutions()[w] = xw.bhp()[w]; - - } - break; - } - } -// const WellType& well_type = wells().type[w]; -// double tot_well_rate = 0.0; -// for (int p = 0; p < np; ++p) { -// tot_well_rate += g[p] * xw.wellRates()[np*w + p]; -// } -// if (well_type == INJECTOR && tot_well_rate == 0) { -// xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; -// xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; -// } - } - } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index b796e80d9..5f5f470b8 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -193,16 +193,16 @@ namespace Opm } // wellSolutions - - if (wells->type[w] == PRODUCER && std::abs(total_well_rates) > 0.0) { - for( int i=0; i 0.0) { + //wellSolutions()[ 0*nw + newIndex ] = prevState.wellSolutions()[0 * nw_old + oldIndex ]; + //if (wells->type[w] == PRODUCER) { + for( int i = 0; i < np; ++i) + { + wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ]; } + //} - - + //} // perfPhaseRates int oldPerf_idx = (*it).second[ 1 ]; @@ -242,6 +242,9 @@ namespace Opm // If the set of controls have changed, this may not be identical // to the last control, but it must be a valid control. currentControls()[ newIndex ] = old_control_index; + WellControls* wc = wells->ctrls[newIndex]; + well_controls_set_current( wc, old_control_index); + } } From 0c9c5a031f5d76f75b0c45f5ddf58fa955322d01 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 19 Aug 2016 10:28:54 +0200 Subject: [PATCH 6/8] Add copy constructure and operator to WellStateFullyImplicitBlackoil --- .../WellStateFullyImplicitBlackoil.hpp | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 5f5f470b8..6de120ba8 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -331,6 +331,26 @@ namespace Opm return res; } + WellStateFullyImplicitBlackoil() = default; + WellStateFullyImplicitBlackoil( const WellStateFullyImplicitBlackoil& rhs ) : + BaseType(rhs), + perfphaserates_( rhs.perfphaserates_ ), + current_controls_( rhs.current_controls_ ), + well_potentials_( rhs.well_potentials_ ), + well_solutions_( rhs.well_solutions_ ) + {} + + WellStateFullyImplicitBlackoil& operator=( const WellStateFullyImplicitBlackoil& rhs ) { + + BaseType::operator =(rhs); + this->perfPhaseRates() = rhs.perfPhaseRates(); + this->currentControls() = rhs.currentControls(); + this->wellPotentials() = rhs.wellPotentials(); + this->wellSolutions() = rhs.wellSolutions(); + return *this; + } + + private: std::vector perfphaserates_; std::vector current_controls_; From cd76816ea59d4bb0a0f56f730c375b26e745eeb8 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 19 Aug 2016 10:30:49 +0200 Subject: [PATCH 7/8] Stay in sync with changes in restartWriter in the eclipseState --- opm/autodiff/FlowMainEbos.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 016af8596..0b2af6a5f 100644 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -66,8 +66,8 @@ namespace Opm // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step) if (Base::param_.has("output_interval")) { const int output_interval = Base::param_.get("output_interval"); - IOConfigPtr ioConfig = Base::eclipse_state_->getIOConfig(); - ioConfig->overrideRestartWriteInterval(static_cast(output_interval)); + eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) ); + } // Possible to force initialization only behavior (NOSIM). From b023cb15a51fc998e273ddbcaa7b0e32d259b7f8 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 19 Aug 2016 13:50:27 +0200 Subject: [PATCH 8/8] Reset wellcontrols from well_state --- opm/autodiff/BlackoilModelEbos.hpp | 2 ++ opm/autodiff/StandardWellsDense.hpp | 9 +++++++++ 2 files changed, 11 insertions(+) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index f72d42bc1..1145a2320 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -283,8 +283,10 @@ namespace Opm { { using namespace Opm::AutoDiffGrid; + //wellModel().assembleAndUpdate(well_state); // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells + wellModel().resetWellControlFromState(well_state); wellModel().updateWellControls(well_state); // Set the primary variables for the wells diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 86ec1034b..71c9052a7 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -86,6 +86,15 @@ namespace Opm { int numPhases() const { return wells().number_of_phases; }; + template + void resetWellControlFromState(WellState xw) { + const int nw = wells_->number_of_wells; + for (int w = 0; w < nw; ++w) { + WellControls* wc = wells_->ctrls[w]; + well_controls_set_current( wc, xw.currentControls()[w]); + } + } + const Wells& wells() const; const Wells* wellsPointer() const;