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

Prepare for extended models.

Browse Source 
Let the code loop over number of components instead of phase
Pass TypeTag as template parameter instead of all the properties.
This commit is contained in:
Tor Harald Sandve 2017-05-03 13:34:15 +02:00
parent 4aa9dc5375
commit 6084721812
10 changed files with 443 additions and 447 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

184
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -118,7 +118,7 @@ namespace Opm {
typedef BlackoilModelParameters ModelParameters;
typedef typename TTAG(EclFlowProblem) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator ;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
@ -129,9 +129,9 @@ namespace Opm {
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef double Scalar;
static const int blocksize = 3;
typedef Dune::FieldVector<Scalar, blocksize > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, blocksize, blocksize > MatrixBlockType;
static const int numEq = BlackoilIndices::numEq;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
@ -162,7 +162,7 @@ namespace Opm {
const ModelParameters& param,
const BlackoilPropsAdFromDeck& fluid,
const DerivedGeology& geo ,
const StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>& well_model,
const StandardWellsDense<TypeTag>& well_model,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output)
: ebosSimulator_(ebosSimulator)
@ -443,7 +443,7 @@ namespace Opm {
{
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int nw = numWells();
return numPhases() * (nc + nw);
return numComponents() * (nc + nw);
}
/// Number of linear iterations used in last call to solveJacobianSystem().
@ -766,7 +766,7 @@ namespace Opm {
template <class CollectiveCommunication>
double convergenceReduction(const CollectiveCommunication& comm,
const long int ncGlobal,
const int np,
const int numComp,
const std::vector< std::vector< Scalar > >& B,
const std::vector< std::vector< Scalar > >& tempV,
const std::vector< std::vector< Scalar > >& R,
@ -777,14 +777,14 @@ namespace Opm {
std::vector< Scalar >& B_avg,
std::vector< Scalar >& maxNormWell )
{
const int nw = residual_well.size() / np;
assert(nw * np == int(residual_well.size()));
const int nw = residual_well.size() / numComp;
assert(nw * numComp == int(residual_well.size()));
// Do the global reductions
B_avg.resize(np);
maxCoeff.resize(np);
R_sum.resize(np);
maxNormWell.resize(np);
B_avg.resize(numComp);
maxCoeff.resize(numComp);
R_sum.resize(numComp);
maxNormWell.resize(numComp);
const std::vector<double>* mask = nullptr;
@ -801,20 +801,20 @@ namespace Opm {
#endif
// computation
for ( int idx = 0; idx < np; ++idx )
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(),
B_avg[compIdx] = std::accumulate( B[ compIdx ].begin(), B[ compIdx ].end(),
0.0 ) / double(ncGlobal);
R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(),
R_sum[compIdx] = std::accumulate( R[ compIdx ].begin(), R[ compIdx ].end(),
0.0 );
maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() ));
maxCoeff[compIdx] = *(std::max_element( tempV[ compIdx ].begin(), tempV[ compIdx ].end() ));
assert(np >= np);
if (idx < np) {
maxNormWell[idx] = 0.0;
assert(numComp >= numComp);
if (compIdx < numComp) {
maxNormWell[compIdx] = 0.0;
for ( int w = 0; w < nw; ++w ) {
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
maxNormWell[compIdx] = std::max(maxNormWell[compIdx], std::abs(residual_well[nw*compIdx + w]));
}
}
}
@ -829,12 +829,12 @@ namespace Opm {
std::vector< Scalar > maxBuffer;
sumBuffer.reserve( B_avg.size() + R_sum.size() + 1 );
maxBuffer.reserve( maxCoeff.size() + maxNormWell.size() );
for( int idx = 0; idx < np; ++idx )
for( int compIdx = 0; compIdx < numComp; ++compIdx )
{
sumBuffer.push_back( B_avg[ idx ] );
sumBuffer.push_back( R_sum[ idx ] );
maxBuffer.push_back( maxCoeff[ idx ] );
maxBuffer.push_back( maxNormWell[ idx ] );
sumBuffer.push_back( B_avg[ compIdx ] );
sumBuffer.push_back( R_sum[ compIdx ] );
maxBuffer.push_back( maxCoeff[ compIdx ] );
maxBuffer.push_back( maxNormWell[ compIdx ] );
}
// Compute total pore volume
@ -847,14 +847,14 @@ namespace Opm {
comm.max( maxBuffer.data(), maxBuffer.size() );
// restore values to local variables
for( int idx = 0, buffIdx = 0; idx < np; ++idx, ++buffIdx )
for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
{
B_avg[ idx ] = sumBuffer[ buffIdx ];
maxCoeff[ idx ] = maxBuffer[ buffIdx ];
B_avg[ compIdx ] = sumBuffer[ buffIdx ];
maxCoeff[ compIdx ] = maxBuffer[ buffIdx ];
++buffIdx;
R_sum[ idx ] = sumBuffer[ buffIdx ];
maxNormWell[ idx ] = maxBuffer[ buffIdx ];
R_sum[ compIdx ] = sumBuffer[ buffIdx ];
maxNormWell[ compIdx ] = maxBuffer[ buffIdx ];
}
// restore global pore volume
@ -880,23 +880,23 @@ namespace Opm {
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = numPhases();
const auto& pv = geo_.poreVolume();
Vector R_sum(np);
Vector B_avg(np);
Vector maxCoeff(np);
Vector maxNormWell(np);
const int numComp = numComponents();
Vector R_sum(numComp);
Vector B_avg(numComp);
Vector maxCoeff(numComp);
Vector maxNormWell(numComp);
// As we will not initialize values in the following arrays
// for the non-interior elements, we have to make sure
// (at least for tempV) that the values there do not influence
// our reduction.
std::vector< Vector > B( np, Vector( nc, 0.0) );
//std::vector< Vector > R( np, Vector( nc, 0.0) ) );
std::vector< Vector > R2( np, Vector( nc, 0.0 ) );
std::vector< Vector > tempV( np, Vector( nc, std::numeric_limits<double>::lowest() ) );
std::vector< Vector > B( numComp, Vector( nc, 0.0) );
//std::vector< Vector > R( numComp, Vector( nc, 0.0) ) );
std::vector< Vector > R2( numComp, Vector( nc, 0.0 ) );
std::vector< Vector > tempV( numComp, Vector( nc, std::numeric_limits<double>::lowest() ) );
const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
@ -915,23 +915,23 @@ namespace Opm {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
for ( int idx = 0; idx < np; ++idx )
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
Vector& R2_idx = R2[ idx ];
Vector& B_idx = B[ idx ];
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(idx);
Vector& R2_idx = R2[ compIdx ];
Vector& B_idx = B[ compIdx ];
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(compIdx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx);
B_idx [cell_idx] = 1.0 / fs.invB(ebosPhaseIdx).value();
R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx];
}
}
for ( int idx = 0; idx < np; ++idx )
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
//tempV.col(idx) = R2.col(idx).abs()/pv;
Vector& tempV_idx = tempV[ idx ];
Vector& R2_idx = R2[ idx ];
//tempV.col(compIdx) = R2.col(compIdx).abs()/pv;
Vector& tempV_idx = tempV[ compIdx ];
Vector& R2_idx = R2[ compIdx ];
for( int cell_idx = 0; cell_idx < nc; ++cell_idx )
{
tempV_idx[ cell_idx ] = std::abs( R2_idx[ cell_idx ] ) / pv[ cell_idx ];
@ -941,32 +941,30 @@ namespace Opm {
Vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
Vector wellResidual = wellModel().residual();
const double pvSum = convergenceReduction(grid_.comm(), global_nc_, np,
const double pvSum = convergenceReduction(grid_.comm(), global_nc_, numComp,
B, tempV, R2, pv_vector, wellResidual,
R_sum, maxCoeff, B_avg, maxNormWell );
Vector CNV(np);
Vector mass_balance_residual(np);
Vector well_flux_residual(np);
Vector CNV(numComp);
Vector mass_balance_residual(numComp);
Vector well_flux_residual(numComp);
bool converged_MB = true;
bool converged_CNV = true;
bool converged_Well = true;
// Finish computation
for ( int idx = 0; idx < np; ++idx )
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
// Well flux convergence is only for fluid phases, not other materials
// in our current implementation.
assert(np >= np);
if (idx < np) {
well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
}
residual_norms.push_back(CNV[idx]);
well_flux_residual[compIdx] = B_avg[compIdx] * maxNormWell[compIdx];
converged_Well = converged_Well && (well_flux_residual[compIdx] < tol_wells);
residual_norms.push_back(CNV[compIdx]);
}
const bool converged = converged_MB && converged_CNV && converged_Well;
@ -977,20 +975,20 @@ namespace Opm {
if (iteration == 0) {
std::string msg = "Iter";
std::vector< std::string > key( np );
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
std::vector< std::string > key( numComp );
for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
key[ phaseIdx ] = std::toupper( phaseName.front() );
}
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
msg += " MB(" + key[ phaseIdx ] + ") ";
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " MB(" + key[ compIdx ] + ") ";
}
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
msg += " CNV(" + key[ phaseIdx ] + ") ";
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " CNV(" + key[ compIdx ] + ") ";
}
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
msg += " W-FLUX(" + key[ phaseIdx ] + ")";
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " W-FLUX(" + key[ compIdx ] + ")";
}
OpmLog::note(msg);
}
@ -998,31 +996,31 @@ namespace Opm {
const std::streamsize oprec = ss.precision(3);
const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
ss << std::setw(4) << iteration;
for (int idx = 0; idx < np; ++idx) {
ss << std::setw(11) << mass_balance_residual[idx];
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << mass_balance_residual[compIdx];
}
for (int idx = 0; idx < np; ++idx) {
ss << std::setw(11) << CNV[idx];
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << CNV[compIdx];
}
for (int idx = 0; idx < np; ++idx) {
ss << std::setw(11) << well_flux_residual[idx];
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << well_flux_residual[compIdx];
}
ss.precision(oprec);
ss.flags(oflags);
OpmLog::note(ss.str());
}
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
for (int phaseIdx = 0; phaseIdx < numPhases(); ++phaseIdx) {
const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
if (std::isnan(mass_balance_residual[phaseIdx])
|| std::isnan(CNV[phaseIdx])
|| (phaseIdx < np && std::isnan(well_flux_residual[phaseIdx]))) {
|| (phaseIdx < numPhases() && std::isnan(well_flux_residual[phaseIdx]))) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName);
}
if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
|| CNV[phaseIdx] > maxResidualAllowed()
|| (phaseIdx < np && well_flux_residual[phaseIdx] > maxResidualAllowed())) {
|| (phaseIdx < numPhases() && well_flux_residual[phaseIdx] > maxResidualAllowed())) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName);
}
}
@ -1037,6 +1035,14 @@ namespace Opm {
return fluid_.numPhases();
}
int numComponents() const
{
if (numPhases() == 2) {
return 2;
}
return FluidSystem::numComponents;
}
/// Wrapper required due to not following generic API
template<class T>
std::vector<std::vector<double> >
@ -1403,13 +1409,8 @@ namespace Opm {
intQuants.pvtRegionIndex()).value();
}
}
const size_t max_num_cells_faillog = 20;
if (failed_cells_pb.size() > 0) {
std::stringstream errlog;
@ -1480,7 +1481,7 @@ namespace Opm {
SimulatorReport failureReport_;
// Well Model
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw> well_model_;
StandardWellsDense<TypeTag> well_model_;
/// \brief Whether we print something to std::cout
bool terminal_output_;
@ -1498,9 +1499,9 @@ namespace Opm {
public:
/// return the StandardWells object
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>&
StandardWellsDense<TypeTag>&
wellModel() { return well_model_; }
const StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>&
const StandardWellsDense<TypeTag>&
wellModel() const { return well_model_; }
/// return the Well struct in the StandardWells
@ -1603,12 +1604,14 @@ namespace Opm {
public:
int ebosCompToFlowPhaseIdx( const int compIdx ) const
{
assert(compIdx < 3);
const int compToPhase[ 3 ] = { Oil, Water, Gas };
return compToPhase[ compIdx ];
}
int flowToEbosPvIdx( const int flowPv ) const
{
assert(flowPv < 3);
const int flowToEbos[ 3 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
@ -1619,6 +1622,7 @@ namespace Opm {
int flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
assert(phaseIdx < 3);
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx };
return phaseToComp[ phaseIdx ];
}
@ -1671,7 +1675,7 @@ namespace Opm {
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(canonicalFlowPhaseIdx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(canonicalFlowPhaseIdx);
const double refDens = FluidSystem::referenceDensity(ebosPhaseIdx, pvtRegionIdx);
for( int pvIdx=0; pvIdx<numFlowPhases; ++pvIdx )
for( int pvIdx = 0; pvIdx < numEq; ++pvIdx )
{
(*col)[ebosCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
(*col)[ebosCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;

2
opm/autodiff/MultisegmentWells_impl.hpp
View File

@ -1040,7 +1040,7 @@ namespace Opm
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_->phaseUsage(),
wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 3. Compute pressure deltas

3
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -57,6 +57,7 @@ public:
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
@ -66,7 +67,7 @@ public:
typedef BlackoilModelEbos Model;
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolver<Model> Solver;
typedef StandardWellsDense<FluidSystem, BlackoilIndices,ElementContext,MaterialLaw> WellModel;
typedef StandardWellsDense<TypeTag> WellModel;
/// Initialise from parameters and objects to observe.

110
opm/autodiff/StandardWellsDense.hpp
View File

@ -2,7 +2,7 @@
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 - 2017 Statoil ASA.
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2016 IRIS AS
Copyright 2016 - 2017 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -67,20 +67,33 @@ enum WellVariablePositions {
/// Class for handling the standard well model.
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
template<typename TypeTag>
class StandardWellsDense {
public:
// --------- Types ---------
typedef WellStateFullyImplicitBlackoilDense WellState;
typedef BlackoilModelParameters ModelParameters;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef double Scalar;
static const int blocksize = 3;
typedef Dune::FieldVector<Scalar, blocksize > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, blocksize, blocksize > MatrixBlockType;
static const int numEq = BlackoilIndices::numEq;
static const int numWellEq = numEq; //number of wellEq is the same as numEq in the model
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<double, /*size=*/blocksize*2> EvalWell;
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
// For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter::
@ -102,27 +115,33 @@ enum WellVariablePositions {
long int global_nc);
template <typename Simulator>
/// The number of components in the model.
int numComponents() const
{
if (phase_usage_.num_phases == 2) {
return 2;
}
return FluidSystem::numComponents;
}
SimulatorReport assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state);
template <typename Simulator>
void assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells);
template <typename Simulator>
void
getMobility(const Simulator& ebosSimulator,
const int perf,
const int cell_idx,
std::vector<EvalWell>& mob) const;
template <typename Simulator>
bool allow_cross_flow(const int w, Simulator& ebosSimulator) const;
bool allow_cross_flow(const int w, const Simulator& ebosSimulator) const;
void localInvert(Mat& istlA) const;
@ -177,8 +196,6 @@ enum WellVariablePositions {
/// Diff to bhp for each well perforation.
const std::vector<double>& wellPerforationPressureDiffs() const;
typedef DenseAd::Evaluation<double, /*size=*/blocksize> Eval;
EvalWell extendEval(Eval in) const;
void setWellVariables(const WellState& xw);
@ -187,12 +204,9 @@ enum WellVariablePositions {
void computeAccumWells();
template<typename FluidState>
void
computeWellFlux(const int& w, const double& Tw, const FluidState& fs, const std::vector<EvalWell>& mob_perfcells_dense,
const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
void computeWellFlux(const int& w, const double& Tw, const IntensiveQuantities& intQuants, const std::vector<EvalWell>& mob_perfcells_dense,
const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
template <typename Simulator>
SimulatorReport solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state);
@ -201,23 +215,18 @@ enum WellVariablePositions {
std::vector<double> residual() const;
template <typename Simulator>
bool getWellConvergence(Simulator& ebosSimulator,
const int iteration) const;
template<typename Simulator>
void
computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw);
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw);
template<typename Simulator>
void
computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf) const;
void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf) const;
void updateWellState(const BVector& dwells,
WellState& well_state) const;
@ -227,12 +236,11 @@ enum WellVariablePositions {
void updateWellControls(WellState& xw) const;
/// 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;
void updateListEconLimited(const Schedule& schedule,
const int current_step,
const Wells* wells_struct,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
void computeWellConnectionDensitesPressures(const WellState& xw,
const std::vector<double>& b_perf,
@ -245,15 +253,12 @@ enum WellVariablePositions {
// Calculating well potentials for each well
// TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
template<typename Simulator>
void
computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
// TODO: some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
template<typename Simulator>
void prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state);
@ -320,14 +325,14 @@ enum WellVariablePositions {
// protected methods
EvalWell getBhp(const int wellIdx) const;
EvalWell getQs(const int wellIdx, const int phaseIdx) const;
EvalWell getQs(const int wellIdx, const int compIdx) const;
EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const;
EvalWell wellVolumeFraction(const int wellIdx, const int compIdx) const;
EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const;
EvalWell wellVolumeFractionScaled(const int wellIdx, const int compIdx) const;
// Q_p / (Q_w + Q_g + Q_o) for three phase cases.
EvalWell wellSurfaceVolumeFraction(const int well_index, const int phase) const;
EvalWell wellSurfaceVolumeFraction(const int well_index, const int compIdx) const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
@ -367,17 +372,14 @@ enum WellVariablePositions {
// TODO: maybe we should provide a light version of computeWellFlux, which does not include the
// calculation of the derivatives
template<typename Simulator>
void
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const int well_index,
std::vector<double>& well_flux) const;
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const int well_index,
std::vector<double>& well_flux) const;
double mostStrictBhpFromBhpLimits(const int well_index) const;
// TODO: maybe it should be improved to be calculate general rates for THP control later
template<typename Simulator>
std::vector<double>
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const int well_index,

522
opm/autodiff/StandardWellsDense_impl.hpp
View File

File diff suppressed because it is too large Load Diff

2
opm/autodiff/StandardWells_impl.hpp
View File

@ -315,7 +315,7 @@ namespace Opm
// Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_->phaseUsage(),
wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];

48
opm/autodiff/WellDensitySegmented.cpp
View File

@ -28,12 +28,10 @@
std::vector<double>
Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& perfComponentRates,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
@ -43,16 +41,17 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
const int numComponents = perfComponentRates.size() / nperf;
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())) {
if (nperf*numComponents != int(surf_dens_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. surf_dens.");
}
if (nperf*np != int(wstate.perfPhaseRates().size())) {
if (nperf*numComponents != int(perfComponentRates.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
}
if (nperf*np != int(b_perf.size())) {
if (nperf*numComponents != int(b_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf.");
}
if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) {
@ -69,20 +68,20 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
// 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);
std::vector<double> q_out_perf(nperf*numComponents);
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) {
for (int component = 0; component < numComponents; ++component) {
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;
q_out_perf[perf*numComponents + component] = 0.0;
} else {
// Set equal to flow from below.
q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
q_out_perf[perf*numComponents + component] = q_out_perf[(perf+1)*numComponents + component];
}
// Subtract outflow through perforation.
q_out_perf[perf*np + phase] -= wstate.perfPhaseRates()[perf*np + phase];
q_out_perf[perf*numComponents + component] -= perfComponentRates[perf*numComponents + component];
}
}
}
@ -93,22 +92,25 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
// 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> mix(numComponents,0.0);
std::vector<double> x(numComponents);
std::vector<double> surf_dens(numComponents);
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);
const double tot_surf_rate = std::accumulate(q_out_perf.begin() + numComponents*perf,
q_out_perf.begin() + numComponents*(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);
for (int component = 0; component < numComponents; ++component) {
mix[component] = std::fabs(q_out_perf[perf*numComponents + component]/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());
for (int phase = 0; phase < np; ++phase) {
mix[phase] = wells.comp_frac[w*np + phase];
}
// intialize 0.0 for comIdx >= np;
}
// Compute volume ratio.
x = mix;
@ -129,11 +131,11 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
x[oilpos] = (mix[oilpos] - mix[gaspos]*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 component = 0; component < numComponents; ++component) {
volrat += x[component] / b_perf[perf*numComponents + component];
}
for (int phase = 0; phase < np; ++phase) {
surf_dens[phase] = surf_dens_perf[perf*np + phase];
for (int component = 0; component < numComponents; ++component) {
surf_dens[component] = surf_dens_perf[perf*numComponents + component];
}
// Compute segment density.

12
opm/autodiff/WellDensitySegmented.hpp
View File

@ -41,17 +41,17 @@ namespace Opm
{
public:
/// Compute well segment densities
/// Notation: N = number of perforations, P = number of phases.
/// Notation: N = number of perforations, C = number of components.
/// \param[in] wells struct with static well info
/// \param[in] wstate dynamic well solution information, only perfRates() is used
/// \param[in] well_rates well rates for actiev components, size NC, P values per perforation
/// \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] b_perf inverse ('little b') formation volume factor, size NC, 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
/// \param[in] surf_dens surface densities for active components, size NC, C values per perforation
static std::vector<double> computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& perfComponentRates,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
@ -80,7 +80,7 @@ namespace Opm
/// Compute pressure deltas.
/// Notation: N = number of perforations, P = number of phases.
/// Notation: N = number of perforations
/// \param[in] wells struct with static well info
/// \param[in] z_perf depth values for each perforation, size N
/// \param[in] dens_perf densities for each perforation, size N (typically computed using computeConnectionDensities)

5
opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
View File

@ -89,7 +89,6 @@ namespace Opm
return;
}
const int np = wells_->number_of_phases;
well_solutions_.clear();
well_solutions_.resize(nw * np, 0.0);
std::vector<double> g = {1.0,1.0,0.01};
@ -130,7 +129,7 @@ namespace Opm
const int waterpos = pu.phase_pos[Water];
const int gaspos = pu.phase_pos[Gas];
assert(np == 3 || (np == 2 && !pu.phase_used[Gas]));
assert(np > 2 || (np == 2 && !pu.phase_used[Gas]));
// assumes the gas fractions are stored after water fractions
if(std::abs(total_rates) > 0) {
if( pu.phase_used[Water] ) {
@ -139,6 +138,8 @@ namespace Opm
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + gaspos] / total_rates ;
}
} else {
if( pu.phase_used[Water] ) {
wellSolutions()[nw + w] = wells_->comp_frac[np*w + waterpos];

2
tests/test_welldensitysegmented.cpp
View File

@ -103,7 +103,7 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
*wells, wellstate, pu,
*wells, pu, rates,
b_perf, rsmax_perf, rvmax_perf, surf_dens);
std::vector<double> dp =
WellDensitySegmented::computeConnectionPressureDelta(