OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 10:40:21 -06:00
Code Issues Packages Projects Releases Wiki Activity

Implement new well model

Browse Source 
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.
This commit is contained in:
Tor Harald Sandve 2016-08-11 13:31:52 +02:00
parent 3e8792b9cc
commit 5cd7468b51
7 changed files with 2684 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -28,7 +28,7 @@
#include <ewoms/common/start.hh>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/StandardWells.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
@ -41,6 +41,7 @@
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
@ -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<ADB> mob_perfcells;
std::vector<ADB> 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<ADB> 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<PhasePresence> 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<int> 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<EigenMatrix> poJac(numPhases + 2);
std::vector<EigenMatrix> poJac(numPhases + 1);
//std::vector<EigenMatrix> TJac(numPhases + 2);
std::vector<std::vector<EigenMatrix>> SJac(numPhases);
std::vector<std::vector<EigenMatrix>> mobJac(numPhases);
std::vector<std::vector<EigenMatrix>> bJac(numPhases);
std::vector<std::vector<EigenMatrix>> pJac(numPhases);
std::vector<EigenMatrix> RsJac(numPhases + 2);
std::vector<EigenMatrix> RvJac(numPhases + 2);
std::vector<EigenMatrix> RsJac(numPhases + 1);
std::vector<EigenMatrix> 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<ADBJacobianMatrix> RsAdbJacs;
std::vector<ADBJacobianMatrix> 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<std::vector<ADBJacobianMatrix>> 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<ADB>& mob_perfcells,
const std::vector<ADB>& 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<V> vars0;
vars0.reserve(2);
vars0.reserve(1);
wellModel().variableWellStateInitials(well_state, vars0);
std::vector<ADB> 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<ADB> 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<M>& Jn = total_residual.derivative();
typedef Eigen::SparseMatrix<double> Sp;
@ -1615,7 +1618,11 @@ namespace Opm {
std::vector<ADB::M> 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<const V>(& well_state.wellSolutions()[0], well_state.wellSolutions().size());
std::vector<ADB::M> old_derivs = state.wellVariables.derivative();
state.wellVariables = ADB::function(std::move(new_well_var), std::move(old_derivs));
//computeWellConnectionPressures(state, well_state);
}
if (!converged) {

2
opm/autodiff/DefaultBlackoilSolutionState.hpp
View File

@ -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<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
};

12
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -390,7 +390,7 @@ namespace Opm
//const int np = residual.material_balance_eq.size();
assert( np == int(residual.material_balance_eq.size()) );
std::vector<ADB> 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;

18
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -25,11 +25,15 @@
#include <opm/autodiff/NonlinearSolver.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
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<Model> 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;

298
opm/autodiff/StandardWellsDense.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_STANDARDWELLSDENSE_HEADER_INCLUDED
#define OPM_STANDARDWELLSDENSE_HEADER_INCLUDED
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <cassert>
#include <tuple>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/DynamicListEconLimited.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
namespace Opm {
/// Class for handling the standard well model.
class StandardWellsDense {
public:
struct WellOps {
explicit WellOps(const Wells* wells);
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
std::vector<int> well_cells; // the set of perforated cells
};
// --------- Types ---------
using ADB = AutoDiffBlock<double>;
typedef DenseAd::Evaluation<double, /*size=*/6> Eval;
//typedef AutoDiffBlock<double> ADB;
using Vector = ADB::V;
using V = ADB::V;
// copied from BlackoilModelBase
// should put to somewhere better
using DataBlock = Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor>;
// --------- Public methods ---------
explicit StandardWellsDense(const Wells* wells_arg);
void init(const BlackoilPropsAdInterface* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* 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 <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
template <class SolutionState>
void computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const;
template <class SolutionState>
void
computeWellFluxDense(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
std::vector<ADB>& 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<Eval>& local, const std::vector<int>& well_cells, const int nc, const std::vector<int>& well_id, const int nw, const int numVars) const;
template <class SolutionState, class WellState>
void updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
template <class WellState>
void updateWellState(const Vector& dwells,
const double dpmaxrel,
WellState& well_state);
template <class WellState>
void updateWellControls(WellState& xw) const;
// TODO: should LinearisedBlackoilResidual also be a template class?
template <class SolutionState>
void addWellFluxEq(const std::vector<ADB>& 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 <class SolutionState, class WellState>
void addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
LinearisedBlackoilResidual& residual);
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
// state0 is non-constant, while it will not be used outside of the function
template <class SolutionState, class WellState>
void
computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state0,
WellState& well_state);
template <class SolutionState, class WellState>
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state,
WellState& well_state) const;
void
variableStateWellIndices(std::vector<int>& indices,
int& next) const;
std::vector<int>
variableWellStateIndices() const;
template <class WellState>
void
variableWellStateInitials(const WellState& xw,
std::vector<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<class WellState>
void
updateListEconLimited(ScheduleConstPtr schedule,
const int current_step,
const Wells* wells,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
template <class SolutionState>
std::vector<ADB> wellVolumeFractions(const SolutionState& state) const;
template <class SolutionState>
void computeAccumWells(const SolutionState& state);
protected:
bool wells_active_;
const Wells* wells_;
const WellOps wops_;
const BlackoilPropsAdInterface* fluid_;
const std::vector<bool>* active_;
const std::vector<PhasePresence>* 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<ADB> F0_;
// protected methods
template <class SolutionState, class WellState>
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf);
template <class WellState>
void computeWellConnectionDensitesPressures(const WellState& xw,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf,
const std::vector<double>& depth_perf,
const double grav);
template <class WellState>
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<bool, bool, int, double>;
enum ConnectionIndex {
INVALIDCONNECTION = -10000
};
template <class WellState>
RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
const WellMapEntryType& map_entry) const;
template <class WellState>
RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
const WellMapEntryType& map_entry) const;
};
} // namespace Opm
#include "StandardWellsDense_impl.hpp"
#endif

2237
opm/autodiff/StandardWellsDense_impl.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

82
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -24,6 +24,7 @@
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <vector>
#include <cassert>
@ -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<double> 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<np; ++i, ++idx, ++oldidx )
{
total_well_rates += prevState.wellRates()[ oldidx ];
}
//if (std::abs(total_well_rates) > 0) {
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
}
// wellSolutions
if (wells->type[w] == PRODUCER && std::abs(total_well_rates) > 0.0) {
for( int i=0; i<np; ++i)
{
wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ];
}
}
// perfPhaseRates
int oldPerf_idx = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ];
@ -149,6 +225,7 @@ namespace Opm
}
}
}
// perfPressures
if( num_perf_old_well == num_perf_this_well )
{
@ -191,6 +268,10 @@ namespace Opm
std::vector<double>& wellPotentials() { return well_potentials_; }
const std::vector<double>& wellPotentials() const { return well_potentials_; }
/// One rate per phase and well connection.
std::vector<double>& wellSolutions() { return well_solutions_; }
const std::vector<double>& wellSolutions() const { return well_solutions_; }
data::Wells report() const override {
data::Wells res = WellState::report();
@ -251,6 +332,7 @@ namespace Opm
std::vector<double> perfphaserates_;
std::vector<int> current_controls_;
std::vector<double> well_potentials_;
std::vector<double> well_solutions_;
};
} // namespace Opm