OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1323 from totto82/refactWell_minimum

Browse Source 
Minimized PR: Make the wellModel self-contained
This commit is contained in:
Atgeirr Flø Rasmussen
2017-11-08 20:52:37 +01:00
committed by GitHub
parent c79dab27d5 a73208f41c
commit a2a0d32e49
16 changed files with 681 additions and 882 deletions
Download Patch File Download Diff File Expand all files Collapse all files

109
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -29,15 +29,8 @@
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
@@ -136,7 +129,7 @@ namespace Opm {
// For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
typedef Opm::FIPData FIPDataType;
@@ -155,7 +148,6 @@ namespace Opm {
BlackoilModelEbos(Simulator& ebosSimulator,
const ModelParameters& param,
BlackoilWellModel<TypeTag>& well_model,
RateConverterType& rate_converter,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output
)
@@ -163,9 +155,6 @@ namespace Opm {
, grid_(ebosSimulator_.gridManager().grid())
, istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
, phaseUsage_(phaseUsageFromDeck(eclState()))
, vfp_properties_(
eclState().getTableManager().getVFPInjTables(),
eclState().getTableManager().getVFPProdTables())
, active_(detail::activePhases(phaseUsage_))
, has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil())
@@ -174,19 +163,12 @@ namespace Opm {
, param_( param )
, well_model_ (well_model)
, terminal_output_ (terminal_output)
, rate_converter_(rate_converter)
, current_relaxation_(1.0)
, dx_old_(AutoDiffGrid::numCells(grid_))
, isBeginReportStep_(false)
{
// Wells are active if they are active wells on at least
// one process.
int wellsActive = localWellsActive() ? 1 : 0;
wellsActive = grid_.comm().max(wellsActive);
wellModel().setWellsActive( wellsActive );
// compute global sum of number of cells
global_nc_ = detail::countGlobalCells(grid_);
wellModel().setVFPProperties(&vfp_properties_);
if (!istlSolver_)
{
OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
@@ -207,13 +189,12 @@ namespace Opm {
const ReservoirState& /*reservoir_state*/,
const WellState& /* well_state */)
{
if ( wellModel().wellCollection()->havingVREPGroups() ) {
updateRateConverter();
}
unsigned numDof = ebosSimulator_.model().numGridDof();
wasSwitched_.resize(numDof);
std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
wellModel().beginTimeStep();
}
@@ -231,7 +212,7 @@ namespace Opm {
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver,
ReservoirState& /*reservoir_state*/,
WellState& well_state)
WellState& /*well_state*/)
{
SimulatorReport report;
failureReport_ = SimulatorReport();
@@ -249,7 +230,7 @@ namespace Opm {
report.total_linearizations = 1;
try {
report += assemble(timer, iteration, well_state);
report += assemble(timer, iteration);
report.assemble_time += perfTimer.stop();
}
catch (...) {
@@ -266,8 +247,8 @@ namespace Opm {
report.converged = getConvergence(timer, iteration,residual_norms) && iteration > nonlinear_solver.minIter();
// checking whether the group targets are converged
if (wellModel().wellCollection()->groupControlActive()) {
report.converged = report.converged && wellModel().wellCollection()->groupTargetConverged(well_state.wellRates());
if (wellModel().wellCollection().groupControlActive()) {
report.converged = report.converged && wellModel().wellCollection().groupTargetConverged(wellModel().wellState().wellRates());
}
report.update_time += perfTimer.stop();
@@ -282,7 +263,6 @@ namespace Opm {
// Compute the nonlinear update.
const int nc = AutoDiffGrid::numCells(grid_);
const int nw = numWells();
BVector x(nc);
try {
@@ -304,11 +284,7 @@ namespace Opm {
// handling well state update before oscillation treatment is a decision based
// on observation to avoid some big performance degeneration under some circumstances.
// there is no theorectical explanation which way is better for sure.
if( nw > 0 )
{
wellModel().recoverWellSolutionAndUpdateWellState(x, well_state);
}
wellModel().recoverWellSolutionAndUpdateWellState(x);
if (param_.use_update_stabilization_) {
// Stabilize the nonlinear update.
@@ -356,6 +332,9 @@ namespace Opm {
DUNE_UNUSED_PARAMETER(timer);
DUNE_UNUSED_PARAMETER(reservoir_state);
DUNE_UNUSED_PARAMETER(well_state);
wellModel().timeStepSucceeded();
}
/// Assemble the residual and Jacobian of the nonlinear system.
@@ -363,18 +342,8 @@ namespace Opm {
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
SimulatorReport assemble(const SimulatorTimerInterface& timer,
const int iterationIdx,
WellState& well_state)
const int iterationIdx)
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// when having VREP group control, update the rate converter based on reservoir state
if ( wellModel().wellCollection()->havingVREPGroups() ) {
updateRateConverter();
}
// -------- Mass balance equations --------
assembleMassBalanceEq(timer, iterationIdx);
@@ -383,7 +352,9 @@ namespace Opm {
try
{
report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
// assembles the well equations and applies the wells to
// the reservoir equations as a source term.
wellModel().assemble(iterationIdx, dt);
// apply well residual to the residual.
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
@@ -394,7 +365,7 @@ namespace Opm {
OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations");
}
return report;
return wellModel().lastReport();
}
// compute the "relative" change of the solution between time steps
@@ -476,14 +447,6 @@ namespace Opm {
}
/// The size (number of unknowns) of the nonlinear system of equations.
int sizeNonLinear() const
{
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int nw = numWells();
return numComponents() * (nc + nw);
}
/// Number of linear iterations used in last call to solveJacobianSystem().
int linearIterationsLastSolve() const
{
@@ -497,6 +460,18 @@ namespace Opm {
const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
// J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
// with the reservoir and D is the wells itself.
// The full system is reduced to a number of cells X number of cells system via Schur complement
// A -= B^T D^-1 C
// Instead of modifying A, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter.
// The residual is modified similarly.
// r = [r, r_well], where r is the residual and r_well the well residual.
// r -= B^T * D^-1 r_well
// Note: Currently the residual is modified in the assembleMassBalanceEq call.
// It conceptually belongs here, but moving it changes the solution path for some
// unknown reasons.
// set initial guess
x = 0.0;
@@ -504,14 +479,14 @@ namespace Opm {
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
Operator opA(ebosJac, wellModel(), istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, well_model_);
Operator opA(ebosJac, wellModel());
istlSolver().solve( opA, x, ebosResid );
}
}
@@ -916,7 +891,7 @@ namespace Opm {
residual_norms.push_back(CNV[compIdx]);
}
const bool converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
const bool converged_Well = wellModel().getWellConvergence(B_avg);
bool converged = converged_MB && converged_Well;
@@ -1486,7 +1461,6 @@ namespace Opm {
const Grid& grid_;
const ISTLSolverType* istlSolver_;
const PhaseUsage phaseUsage_;
VFPProperties vfp_properties_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active phases. Maps active -> canonical phase indices.
@@ -1507,9 +1481,6 @@ namespace Opm {
/// \brief The number of cells of the global grid.
long int global_nc_;
// rate converter between the surface volume rates and reservoir voidage rates
RateConverterType& rate_converter_;
std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_;
BVector dx_old_;
@@ -1523,15 +1494,6 @@ namespace Opm {
const BlackoilWellModel<TypeTag>&
wellModel() const { return well_model_; }
int numWells() const { return well_model_.numWells(); }
/// return true if wells are available on this process
bool localWellsActive() const { return well_model_.localWellsActive(); }
public:
int flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
const auto& pu = phaseUsage_;
@@ -1561,15 +1523,6 @@ namespace Opm {
return phaseIdx;
}
private:
void updateRateConverter()
{
rate_converter_.defineState<ElementContext>(ebosSimulator_);
}
public:
void beginReportStep()
{
isBeginReportStep_ = true;

3
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@@ -116,6 +116,8 @@ namespace Opm {
using Base::numPhases;
using Base::numMaterials;
using Base::materialName;
using Base::wellModel;
protected:
// --------- Data members ---------
@@ -143,7 +145,6 @@ namespace Opm {
using Base::vfp_properties_;
using Base::well_model_;
using Base::wellModel;
// using Base::wells;
using Base::wellsActive;
using Base::updatePrimalVariableFromState;

196
opm/autodiff/BlackoilWellModel.hpp
View File

@@ -50,6 +50,7 @@
#include <opm/autodiff/WellInterface.hpp>
#include <opm/autodiff/StandardWell.hpp>
#include <opm/autodiff/MultisegmentWell.hpp>
#include<opm/autodiff/SimFIBODetails.hpp>
#include<dune/common/fmatrix.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/matrixmatrix.hh>
@@ -69,6 +70,10 @@ namespace Opm {
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
static const int Water = WellInterface<TypeTag>::Water;
static const int Oil = WellInterface<TypeTag>::Oil;
static const int Gas = WellInterface<TypeTag>::Gas;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
@@ -88,32 +93,16 @@ namespace Opm {
// For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
// --------- Public methods ---------
BlackoilWellModel(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
BlackoilWellModel(Simulator& ebosSimulator,
const ModelParameters& param,
const RateConverterType& rate_converter,
const bool terminal_output,
const int current_index,
const std::vector<int>& pvt_region_idx);
const bool terminal_output);
void init(const PhaseUsage phase_usage_arg,
const std::vector<bool>& active_arg,
const double gravity_arg,
const std::vector<double>& depth_arg,
long int global_nc,
const Grid& grid);
void setVFPProperties(const VFPProperties* vfp_properties_arg);
SimulatorReport assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state);
// compute the well fluxes and assemble them in to the reservoir equations as source terms
// and in the well equations.
void assemble(const int iterationIdx,
const double dt);
// substract Binv(D)rw from r;
void apply( BVector& r) const;
@@ -126,44 +115,46 @@ namespace Opm {
// using the solution x to recover the solution xw for wells and applying
// xw to update Well State
void recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const;
void recoverWellSolutionAndUpdateWellState(const BVector& x);
int numWells() const;
// Check if well equations is converged.
bool getWellConvergence(const std::vector<Scalar>& B_avg) const;
/// return true if wells are available in the reservoir
bool wellsActive() const;
// return all the wells.
const WellCollection& wellCollection() const;
// return non const reference to all the wells.
WellCollection& wellCollection();
void setWellsActive(const bool wells_active);
// return the internal well state, ignore the passed one.
// Used by the legacy code to make it compatible with the legacy well models.
const WellState& wellState(const WellState& well_state OPM_UNUSED) const;
/// return true if wells are available on this process
bool localWellsActive() const;
// return the internal well state
const WellState& wellState() const;
bool getWellConvergence(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const;
// only use this for restart.
void setRestartWellState(const WellState& well_state);
/// upate the dynamic lists related to economic limits
void updateListEconLimited(const Schedule& schedule,
const int current_step,
const Wells* wells_struct,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
// called at the beginning of a time step
void beginTimeStep();
// called at the end of a time step
void timeStepSucceeded();
WellCollection* wellCollection() const;
// called at the beginning of a report step
void beginReportStep(const int time_step);
// called at the end of a report step
void endReportStep();
const SimulatorReport& lastReport() const;
protected:
Simulator& ebosSimulator_;
std::unique_ptr<WellsManager> wells_manager_;
std::vector< const Well* > wells_ecl_;
bool wells_active_;
const Wells* wells_;
const std::vector< const Well* > wells_ecl_;
// the number of wells in this process
// trying not to use things from Wells struct
// TODO: maybe a better name to emphasize it is local?
const int number_of_wells_;
const int number_of_phases_;
const ModelParameters& param_;
using WellInterfacePtr = std::unique_ptr<WellInterface<TypeTag> >;
// a vector of all the wells.
@@ -175,59 +166,60 @@ namespace Opm {
using ConvergenceReport = typename WellInterface<TypeTag>::ConvergenceReport;
// create the well container
static std::vector<WellInterfacePtr > createWellContainer(const Wells* wells,
const std::vector<const Well*>& wells_ecl,
const bool use_multisegment_well,
const int time_step,
const ModelParameters& param);
std::vector<WellInterfacePtr > createWellContainer(const int time_step) const;
// Well collection is used to enforce the group control
WellCollection* well_collection_;
WellState well_state_;
WellState previous_well_state_;
const ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;
int current_timeIdx_;
std::vector<int> pvt_region_idx_;
PhaseUsage phase_usage_;
std::vector<bool> active_;
const RateConverterType& rate_converter_;
const std::vector<int>& pvt_region_idx_;
size_t global_nc_;
// the number of the cells in the local grid
int number_of_cells_;
size_t number_of_cells_;
double gravity_;
std::vector<double> depth_;
long int global_nc_;
DynamicListEconLimited dynamic_list_econ_limited_;
std::unique_ptr<RateConverterType> rateConverter_;
std::unique_ptr<VFPProperties> vfp_properties_;
SimulatorReport last_report_;
// used to better efficiency of calcuation
mutable BVector scaleAddRes_;
void updateWellControls(WellState& xw) const;
const Wells* wells() const { return wells_manager_->c_wells(); }
void updateGroupControls(WellState& well_state) const;
const Schedule& schedule() const
{ return ebosSimulator_.gridManager().schedule(); }
void updateWellControls();
void updateGroupControls();
// setting the well_solutions_ based on well_state.
void updatePrimaryVariables(const WellState& well_state) const;
void updatePrimaryVariables();
void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
void computeRepRadiusPerfLength(const Grid& grid);
void computeAverageFormationFactor(const Simulator& ebosSimulator,
std::vector<double>& B_avg) const;
void computeAverageFormationFactor(std::vector<double>& B_avg) const;
void applyVREPGroupControl(WellState& well_state) const;
void applyVREPGroupControl();
void computeWellVoidageRates(const WellState& well_state,
std::vector<double>& well_voidage_rates,
void computeWellVoidageRates(std::vector<double>& well_voidage_rates,
std::vector<double>& voidage_conversion_coeffs) const;
// Calculating well potentials for each well
// TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
void computeWellPotentials(std::vector<double>& well_potentials);
const std::vector<double>& wellPerfEfficiencyFactors() const;
@@ -238,47 +230,51 @@ namespace Opm {
// twice at the beginning of the time step
/// Calculating the explict quantities used in the well calculation. By explicit, we mean they are cacluated
/// at the beginning of the time step and no derivatives are included in these quantities
void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& xw) const;
void calculateExplicitQuantities() const;
SimulatorReport solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state) const;
SimulatorReport solveWellEq(const double dt);
void initPrimaryVariablesEvaluation() const;
// The number of components in the model.
int numComponents() const
{
if (numPhases() == 2) {
return 2;
}
int numComp = FluidSystem::numComponents;
if (has_solvent_) {
numComp ++;
}
int numComponents() const;
return numComp;
}
int numWells() const;
int numPhases() const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
void resetWellControlFromState(const WellState& xw) const;
void resetWellControlFromState() const;
void assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells) const;
void assembleWellEq(const double dt,
bool only_wells);
// some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
void prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state) const;
void prepareTimeStep();
void prepareGroupControl(const Simulator& ebos_simulator,
WellState& well_state) const;
void prepareGroupControl();
void computeRESV(const std::size_t step);
void extractLegacyCellPvtRegionIndex_();
void extractLegacyDepth_();
/// 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;
/// upate the dynamic lists related to economic limits
void updateListEconLimited(DynamicListEconLimited& list_econ_limited) const;
void updatePerforationIntensiveQuantities();
};

738
opm/autodiff/BlackoilWellModel_impl.hpp
View File

File diff suppressed because it is too large Load Diff

6
opm/autodiff/FlowMainEbos.hpp
View File

@@ -728,8 +728,7 @@ namespace Opm
FluidSystem::enableDissolvedGas(),
FluidSystem::enableVaporizedOil(),
eclState(),
*output_writer_,
defunctWellNames()));
*output_writer_));
}
private:
@@ -806,9 +805,6 @@ namespace Opm
Scalar gravity() const
{ return ebosProblem().gravity()[2]; }
std::unordered_set<std::string> defunctWellNames() const
{ return ebosSimulator_->gridManager().defunctWellNames(); }
data::Solution computeLegacySimProps_()
{
const int* dims = UgGridHelpers::cartDims(grid());

3
opm/autodiff/MultisegmentWell.hpp
View File

@@ -47,6 +47,9 @@ namespace Opm
using Base::has_solvent;
using Base::has_polymer;
using Base::Water;
using Base::Oil;
using Base::Gas;
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.

4
opm/autodiff/MultisegmentWells.hpp
View File

@@ -237,6 +237,10 @@ namespace Opm {
const Vector& wellPerfEfficiencyFactors() const;
// just return the passed well state
template<class WellState>
const WellState& wellState(const WellState& well_state) const { return well_state; }
protected:
// TODO: probably a wells_active_ will be required here.

3
opm/autodiff/RateConverter.hpp
View File

@@ -21,13 +21,12 @@
#ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
#define OPM_RATECONVERTER_HPP_HEADER_INCLUDED
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>
#include <cmath>
#include <memory>

341
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -99,8 +99,7 @@ public:
const bool has_disgas,
const bool has_vapoil,
const EclipseState& /* eclState */,
OutputWriter& output_writer,
const std::unordered_set<std::string>& defunct_well_names)
OutputWriter& output_writer)
: ebosSimulator_(ebosSimulator),
param_(param),
model_param_(param),
@@ -111,8 +110,6 @@ public:
has_vapoil_(has_vapoil),
terminal_output_(param.getDefault("output_terminal", true)),
output_writer_(output_writer),
rateConverter_(createRateConverter_()),
defunct_well_names_( defunct_well_names ),
is_parallel_run_( false )
{
#if HAVE_MPI
@@ -141,7 +138,6 @@ public:
ExtraData extra;
failureReport_ = SimulatorReport();
extractLegacyDepth_();
// communicate the initial solution to ebos
if (timer.initialStep()) {
@@ -214,7 +210,6 @@ public:
ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
}
DynamicListEconLimited dynamic_list_econ_limited;
SimulatorReport report;
SimulatorReport stepReport;
@@ -231,6 +226,13 @@ public:
std::vector<std::vector<double>> originalFluidInPlace;
std::vector<double> originalFluidInPlaceTotals;
WellModel well_model(ebosSimulator_, model_param_, terminal_output_);
if (output_writer_.isRestart()) {
well_model.setRestartWellState(prev_well_state); // Neccessary for perfect restarts
}
WellState wellStateDummy; //not used. Only passed to make the old interfaces happy
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@@ -242,77 +244,10 @@ public:
OpmLog::note(ss.str());
}
// Create wells and well state.
WellsManager wells_manager(eclState(),
schedule(),
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(grid()),
Opm::UgGridHelpers::globalCell(grid()),
Opm::UgGridHelpers::cartDims(grid()),
Opm::UgGridHelpers::dimensions(grid()),
Opm::UgGridHelpers::cell2Faces(grid()),
Opm::UgGridHelpers::beginFaceCentroids(grid()),
dynamic_list_econ_limited,
is_parallel_run_,
defunct_well_names_ );
const Wells* wells = wells_manager.c_wells();
WellState well_state;
// The well state initialize bhp with the cell pressure in the top cell.
// We must therefore provide it with updated cell pressures
size_t nc = Opm::UgGridHelpers::numCells(grid());
std::vector<double> cellPressures(nc, 0.0);
const auto& gridView = ebosSimulator_.gridManager().gridView();
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
cellPressures[cellIdx] = p;
}
well_state.init(wells, cellPressures, prev_well_state, phaseUsage_);
// give the polymer and surfactant simulators the chance to do their stuff
handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
// Compute reservoir volumes for RESV controls.
computeRESV(timer.currentStepNum(), wells, well_state);
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
const auto& wells_ecl = schedule().getWells(timer.currentStepNum());
extractLegacyCellPvtRegionIndex_();
WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_,
rateConverter_, terminal_output_, timer.currentStepNum(), legacyCellPvtRegionIdx_);
// handling MS well related
if (model_param_.use_multisegment_well_) { // if we use MultisegmentWell model
for (const auto& well : wells_ecl) {
// TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well
// while keep the well shut. More accurately, we should check if the well exisits in the Wells
// structure here
if (well->isMultiSegment(timer.currentStepNum()) ) { // there is one well is MS well
well_state.initWellStateMSWell(wells, wells_ecl, timer.currentStepNum(), phaseUsage_, prev_well_state);
break;
}
}
}
well_model.beginReportStep(timer.currentStepNum());
auto solver = createSolver(well_model);
@@ -344,7 +279,7 @@ public:
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
output_writer_.writeTimeStep( timer, state, well_model.wellState(), solver->model() );
report.output_write_time += perfTimer.stop();
}
@@ -373,14 +308,14 @@ public:
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
stepReport = adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
stepReport = adaptiveTimeStepping->step( timer, *solver, state, wellStateDummy, event, output_writer_,
output_writer_.requireFIPNUM() ? &fipnum : nullptr );
report += stepReport;
failureReport_ += adaptiveTimeStepping->failureReport();
}
else {
// solve for complete report step
stepReport = solver->step(timer, state, well_state);
stepReport = solver->step(timer, state, wellStateDummy);
report += stepReport;
failureReport_ += solver->failureReport();
@@ -400,6 +335,7 @@ public:
}
solver->model().endReportStep();
well_model.endReportStep();
// take time that was used to solve system for this reportStep
solver_timer.stop();
@@ -450,13 +386,9 @@ public:
Dune::Timer perfTimer;
perfTimer.start();
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
output_writer_.writeTimeStep( timer, state, well_state, solver->model(), false, nextstep, report);
output_writer_.writeTimeStep( timer, state, well_model.wellState(), solver->model(), false, nextstep, report);
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;
updateListEconLimited(solver, schedule(), timer.currentStepNum(), wells,
well_state, dynamic_list_econ_limited);
}
// Stop timer and create timing report
@@ -474,222 +406,18 @@ public:
{ return ebosSimulator_.gridManager().grid(); }
protected:
void handleAdditionalWellInflow(SimulatorTimer& /*timer*/,
WellsManager& /* wells_manager */,
WellState& /* well_state */,
const Wells* /* wells */)
{
}
std::unique_ptr<Solver> createSolver(WellModel& well_model)
{
const auto& gridView = ebosSimulator_.gridView();
const PhaseUsage& phaseUsage = phaseUsage_;
const std::vector<bool> activePhases = detail::activePhases(phaseUsage);
const double gravity = ebosSimulator_.problem().gravity()[2];
// calculate the number of elements of the compressed sequential grid. this needs
// to be done in two steps because the dune communicator expects a reference as
// argument for sum()
int globalNumCells = gridView.size(/*codim=*/0);
globalNumCells = gridView.comm().sum(globalNumCells);
well_model.init(phaseUsage,
activePhases,
gravity,
legacyDepth_,
globalNumCells,
grid());
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
model_param_,
well_model,
rateConverter_,
solver_,
terminal_output_));
return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
}
void computeRESV(const std::size_t step,
const Wells* wells,
WellState& xw)
{
typedef SimFIBODetails::WellMap WellMap;
const auto w_ecl = schedule().getWells(step);
const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
const std::size_t number_resv_wells = resv_wells.size();
std::size_t global_number_resv_wells = number_resv_wells;
#if HAVE_MPI
if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
const auto& info =
boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
if ( global_number_resv_wells )
{
// At least one process has resv wells. Therefore rate converter needs
// to calculate averages over regions that might cross process
// borders. This needs to be done by all processes and therefore
// outside of the next if statement.
rateConverter_.template defineState<ElementContext>(ebosSimulator_);
}
}
else
#endif
{
if ( global_number_resv_wells )
{
rateConverter_.template defineState<ElementContext>(ebosSimulator_);
}
}
if (! resv_wells.empty()) {
const PhaseUsage& pu = phaseUsage_;
const std::vector<double>::size_type np = phaseUsage_.num_phases;
std::vector<double> distr (np);
std::vector<double> hrates(np);
std::vector<double> prates(np);
for (std::vector<int>::const_iterator
rp = resv_wells.begin(), e = resv_wells.end();
rp != e; ++rp)
{
WellControls* ctrl = wells->ctrls[*rp];
const bool is_producer = wells->type[*rp] == PRODUCER;
const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
const auto& eclProblem = ebosSimulator_.problem();
const int pvtreg = eclProblem.pvtRegionIndex(well_cell_top);
// RESV control mode, all wells
{
const int rctrl = SimFIBODetails::resv_control(ctrl);
if (0 <= rctrl) {
const std::vector<double>::size_type off = (*rp) * np;
if (is_producer) {
// Convert to positive rates to avoid issues
// in coefficient calculations.
std::transform(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin(), std::negate<double>());
} else {
std::copy(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin());
}
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(fipreg, pvtreg, distr);
well_controls_iset_distr(ctrl, rctrl, & distr[0]);
}
}
// RESV control, WCONHIST wells. A bit of duplicate
// work, regrettably.
if (is_producer && wells->name[*rp] != 0) {
WellMap::const_iterator i = wmap.find(wells->name[*rp]);
if (i != wmap.end()) {
const auto* wp = i->second;
const WellProductionProperties& p =
wp->getProductionProperties(step);
if (! p.predictionMode) {
// History matching (WCONHIST/RESV)
SimFIBODetails::historyRates(pu, p, hrates);
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(fipreg, pvtreg, distr);
// WCONHIST/RESV target is sum of all
// observed phase rates translated to
// reservoir conditions. Recall sign
// convention: Negative for producers.
const double target =
- std::inner_product(distr.begin(), distr.end(),
hrates.begin(), 0.0);
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
static const double invalid_alq = -std::numeric_limits<double>::max();
static const int invalid_vfp = -std::numeric_limits<int>::max();
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
invalid_alq, invalid_vfp,
& distr[0], ctrl);
// For WCONHIST the BHP limit is set to 1 atm.
// or a value specified using WELTARG
double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
const int ok_bhp =
well_controls_add_new(BHP, bhp_limit,
invalid_alq, invalid_vfp,
NULL, ctrl);
if (ok_resv != 0 && ok_bhp != 0) {
xw.currentControls()[*rp] = 0;
well_controls_set_current(ctrl, 0);
}
}
}
}
}
}
if( wells )
{
for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
WellControls* ctrl = wells->ctrls[w];
const bool is_producer = wells->type[w] == PRODUCER;
if (!is_producer && wells->name[w] != 0) {
WellMap::const_iterator i = wmap.find(wells->name[w]);
if (i != wmap.end()) {
const auto* wp = i->second;
const WellInjectionProperties& injector = wp->getInjectionProperties(step);
if (!injector.predictionMode) {
//History matching WCONINJEH
static const double invalid_alq = -std::numeric_limits<double>::max();
static const int invalid_vfp = -std::numeric_limits<int>::max();
// For WCONINJEH the BHP limit is set to a large number
// or a value specified using WELTARG
double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
const int ok_bhp =
well_controls_add_new(BHP, bhp_limit,
invalid_alq, invalid_vfp,
NULL, ctrl);
if (!ok_bhp) {
OPM_THROW(std::runtime_error, "Failed to add well control.");
}
}
}
}
}
}
}
void updateListEconLimited(const std::unique_ptr<Solver>& solver,
const Schedule& schedule,
const int current_step,
const Wells* wells,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const
{
solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
well_state, list_econ_limited);
}
void FIPUnitConvert(const UnitSystem& units,
std::vector<std::vector<double>>& fip)
{
@@ -845,23 +573,6 @@ protected:
const Schedule& schedule() const
{ return ebosSimulator_.gridManager().schedule(); }
const SummaryConfig& summaryConfig() const
{ return ebosSimulator_.gridManager().summaryConfig( ); }
void extractLegacyCellPvtRegionIndex_()
{
const auto& grid = ebosSimulator_.gridManager().grid();
const auto& eclProblem = ebosSimulator_.problem();
const unsigned numCells = grid.size(/*codim=*/0);
legacyCellPvtRegionIdx_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
legacyCellPvtRegionIdx_[cellIdx] =
eclProblem.pvtRegionIndex(cellIdx);
}
}
void initHysteresisParams(ReservoirState& state) {
const int num_cells = Opm::UgGridHelpers::numCells(grid());
@@ -895,17 +606,6 @@ protected:
}
}
void extractLegacyDepth_()
{
const auto& grid = ebosSimulator_.gridManager().grid();
const unsigned numCells = grid.size(/*codim=*/0);
legacyDepth_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
legacyDepth_[cellIdx] =
grid.cellCenterDepth(cellIdx);
}
}
// Used to convert initial Reservoirstate to primary variables in the SolutionVector
void convertInput( const int iterationIdx,
@@ -1014,18 +714,10 @@ protected:
}
}
RateConverterType createRateConverter_() {
RateConverterType rate_converter(phaseUsage_,
std::vector<int>(AutoDiffGrid::numCells(grid()), 0)); // FIP = 0
return rate_converter;
}
// Data.
Simulator& ebosSimulator_;
std::vector<int> legacyCellPvtRegionIdx_;
std::vector<double> legacyDepth_;
typedef typename Solver::SolverParameters SolverParameters;
SimulatorReport failureReport_;
@@ -1043,11 +735,6 @@ protected:
bool terminal_output_;
// output_writer
OutputWriter& output_writer_;
RateConverterType rateConverter_;
// The names of wells that should be defunct
// (e.g. in a parallel run when they are handeled by
// a different process)
std::unordered_set<std::string> defunct_well_names_;
// Whether this a parallel simulation or not
bool is_parallel_run_;

3
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@@ -1045,8 +1045,7 @@ namespace Opm
miscSummaryData["TCPU"] = totalSolverTime;
}
}
writeTimeStepWithCellProperties(timer, localState, localCellData, localWellState, miscSummaryData, extraRestartData, substep);
writeTimeStepWithCellProperties(timer, localState, localCellData, physicalModel.wellModel().wellState(localWellState), miscSummaryData, extraRestartData, substep);
}
}
#endif

3
opm/autodiff/StandardWell.hpp
View File

@@ -67,6 +67,9 @@ namespace Opm
using Base::has_solvent;
using Base::has_polymer;
using Base::name;
using Base::Water;
using Base::Oil;
using Base::Gas;
// TODO: with flow_ebos，for a 2P deck, // TODO: for the 2p deck, numEq will be 3, a dummy phase is already added from the reservoir side.
// it will cause problem here without processing the dummy phase.

5
opm/autodiff/StandardWells.hpp
View File

@@ -199,6 +199,11 @@ namespace Opm {
const Vector& wellPerfEfficiencyFactors() const;
// just return the passed well state
template<class WellState>
const WellState& wellState(const WellState& well_state) const { return well_state; }
protected:
bool wells_active_;
const Wells* wells_;

120
opm/autodiff/WellDensitySegmented.cpp
View File

@@ -20,7 +20,6 @@
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/core/wells.h>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <numeric>
@@ -148,125 +147,6 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
std::vector<double>
Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoilSolvent& wstate,
const PhaseUsage& phase_usage,
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)
{
// Verify that we have consistent input.
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
if (wells.number_of_phases != phase_usage.num_phases) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. phase_usage.");
}
if (nperf*np != int(surf_dens_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. surf_dens.");
}
if (nperf*np != int(wstate.perfPhaseRates().size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
}
if (nperf*np != int(b_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf.");
}
if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) {
// Need both oil and gas phases.
if (!phase_usage.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::logic_error, "Oil phase inactive, but non-empty rsmax_perf or rvmax_perf.");
}
if (!phase_usage.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::logic_error, "Gas phase inactive, but non-empty rsmax_perf or rvmax_perf.");
}
}
// 1. Compute the flow (in surface volume units for each
// component) exiting up the wellbore from each perforation,
// taking into account flow from lower in the well, and
// in/out-flow at each perforation.
std::vector<double> q_out_perf(nperf*np);
for (int w = 0; w < nw; ++w) {
// Iterate over well perforations from bottom to top.
for (int perf = wells.well_connpos[w+1] - 1; perf >= wells.well_connpos[w]; --perf) {
for (int phase = 0; phase < np; ++phase) {
if (perf == wells.well_connpos[w+1] - 1) {
// This is the bottom perforation. No flow from below.
q_out_perf[perf*np + phase] = 0.0;
} else {
// Set equal to flow from below.
q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
}
// Subtract outflow through perforation.
q_out_perf[perf*np + phase] -= wstate.perfPhaseRates()[perf*np + phase];
}
}
}
// 2. Compute the component mix at each perforation as the
// absolute values of the surface rates divided by their sum.
// Then compute volume ratios (formation factors) for each perforation.
// Finally compute densities for the segments associated with each perforation.
const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
std::vector<double> mix(np);
std::vector<double> x(np);
std::vector<double> surf_dens(np);
std::vector<double> dens(nperf);
for (int w = 0; w < nw; ++w) {
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
// Find component mix.
const double tot_surf_rate = std::accumulate(q_out_perf.begin() + np*perf,
q_out_perf.begin() + np*(perf+1), 0.0);
if (tot_surf_rate != 0.0) {
for (int phase = 0; phase < np; ++phase) {
mix[phase] = std::fabs(q_out_perf[perf*np + phase]/tot_surf_rate);
}
} else {
// No flow => use well specified fractions for mix.
std::copy(wells.comp_frac + w*np, wells.comp_frac + (w+1)*np, mix.begin());
}
// Compute volume ratio.
x = mix;
double rs = 0.0;
double rv = 0.0;
if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
}
const double gas_without_solvent = mix[gaspos] * (1.0 - wstate.solventFraction()[w]);
if (!rvmax_perf.empty() && gas_without_solvent > 0.0) {
rv = std::min(mix[oilpos]/gas_without_solvent, rvmax_perf[perf]);
}
if (rs != 0.0) {
// Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
}
if (rv != 0.0) {
// Subtract oil in gas from oil mixture
x[oilpos] = (mix[oilpos] - gas_without_solvent*rv)/(1.0 - rs*rv);;
}
double volrat = 0.0;
for (int phase = 0; phase < np; ++phase) {
volrat += x[phase] / b_perf[perf*np + phase];
}
for (int phase = 0; phase < np; ++phase) {
surf_dens[phase] = surf_dens_perf[perf*np + phase];
}
// Compute segment density.
dens[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
}
}
return dens;
}
std::vector<double>
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const std::vector<double>& z_perf,

19
opm/autodiff/WellDensitySegmented.hpp
View File

@@ -59,25 +59,6 @@ namespace Opm
/// Compute well segment densities for solvent model
/// Notation: N = number of perforations, P = number of phases.
/// \param[in] wells struct with static well info
/// \param[in] wstate dynamic well solution information, perfRates() and solventFraction() is used
/// \param[in] phase_usage specifies which phases are active and not
/// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation
/// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N
/// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N
/// \param[in] surf_dens surface densities for active components, size NP, P values per perforation
static std::vector<double> computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoilSolvent& wstate,
const PhaseUsage& phase_usage,
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);
/// Compute pressure deltas.
/// Notation: N = number of perforations

5
opm/autodiff/WellInterface.hpp
View File

@@ -64,6 +64,11 @@ namespace Opm
using WellState = WellStateFullyImplicitBlackoil;
typedef BlackoilModelParameters ModelParameters;
static const int Water = BlackoilPhases::Aqua;
static const int Oil = BlackoilPhases::Liquid;
static const int Gas = BlackoilPhases::Vapour;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;

5
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@@ -21,7 +21,6 @@
#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
@@ -375,8 +374,8 @@ namespace Opm
const int start_perf = wells->well_connpos[w];
const int start_perf_next_well = wells->well_connpos[w + 1];
assert(nperf == (start_perf_next_well - start_perf)); // make sure the information from wells_ecl consistent with wells
if (pu.phase_used[Gas]) {
const int gaspos = pu.phase_pos[Gas];
if (pu.phase_used[BlackoilPhases::Vapour]) {
const int gaspos = pu.phase_pos[BlackoilPhases::Vapour];
// scale the phase rates for Gas to avoid too bad initial guess for gas fraction
// it will probably benefit the standard well too, while it needs to be justified
// TODO: to see if this strategy can benefit StandardWell too