From 608472181200f66a6d8435b6ad76970ef3dc53fc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 3 May 2017 13:34:15 +0200 Subject: [PATCH] Prepare for extended models. Let the code loop over number of components instead of phase Pass TypeTag as template parameter instead of all the properties. --- opm/autodiff/BlackoilModelEbos.hpp | 184 +++--- opm/autodiff/MultisegmentWells_impl.hpp | 2 +- .../SimulatorFullyImplicitBlackoilEbos.hpp | 3 +- opm/autodiff/StandardWellsDense.hpp | 110 ++-- opm/autodiff/StandardWellsDense_impl.hpp | 522 +++++++++--------- opm/autodiff/StandardWells_impl.hpp | 2 +- opm/autodiff/WellDensitySegmented.cpp | 48 +- opm/autodiff/WellDensitySegmented.hpp | 12 +- .../WellStateFullyImplicitBlackoilDense.hpp | 5 +- tests/test_welldensitysegmented.cpp | 2 +- 10 files changed, 443 insertions(+), 447 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 31d6a0a45..bcd64719d 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -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 VectorBlockType; - typedef Dune::FieldMatrix MatrixBlockType; + static const int numEq = BlackoilIndices::numEq; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; @@ -162,7 +162,7 @@ namespace Opm { const ModelParameters& param, const BlackoilPropsAdFromDeck& fluid, const DerivedGeology& geo , - const StandardWellsDense& well_model, + const StandardWellsDense& 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 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* 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::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::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 std::vector > @@ -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 well_model_; + StandardWellsDense 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& + StandardWellsDense& wellModel() { return well_model_; } - const StandardWellsDense& + const StandardWellsDense& 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 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 diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 30a39096a..346f860d9 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -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 Solver; - typedef StandardWellsDense WellModel; + typedef StandardWellsDense WellModel; /// Initialise from parameters and objects to observe. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 9d17afcae..5c2906a49 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -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 + template 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 VectorBlockType; - typedef Dune::FieldMatrix 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 VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; - typedef DenseAd::Evaluation EvalWell; + typedef DenseAd::Evaluation EvalWell; + typedef DenseAd::Evaluation 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 + /// 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 void assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, bool only_wells); - template void getMobility(const Simulator& ebosSimulator, const int perf, const int cell_idx, std::vector& mob) const; - template - 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& wellPerforationPressureDiffs() const; - typedef DenseAd::Evaluation Eval; - EvalWell extendEval(Eval in) const; void setWellVariables(const WellState& xw); @@ -187,12 +204,9 @@ enum WellVariablePositions { void computeAccumWells(); - template - void - computeWellFlux(const int& w, const double& Tw, const FluidState& fs, const std::vector& mob_perfcells_dense, - const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const; + void computeWellFlux(const int& w, const double& Tw, const IntensiveQuantities& intQuants, const std::vector& mob_perfcells_dense, + const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const; - template SimulatorReport solveWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state); @@ -201,23 +215,18 @@ enum WellVariablePositions { std::vector residual() const; - template bool getWellConvergence(Simulator& ebosSimulator, const int iteration) const; - template - void - computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw); + void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw); - template - void - computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) const; + void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& 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& 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 - void - computeWellPotentials(const Simulator& ebosSimulator, - const WellState& well_state, - std::vector& well_potentials) const; + void computeWellPotentials(const Simulator& ebosSimulator, + const WellState& well_state, + std::vector& 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 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 - void - computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, - const int well_index, - std::vector& well_flux) const; + void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const EvalWell& bhp, + const int well_index, + std::vector& 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 std::vector computeWellPotentialWithTHP(const Simulator& ebosSimulator, const int well_index, diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 330ff6476..d1af3d5dd 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -3,8 +3,8 @@ namespace Opm { - template - StandardWellsDense:: + template + StandardWellsDense:: StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, const ModelParameters& param, @@ -17,8 +17,8 @@ namespace Opm { , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0) , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) - , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) - , F0_(wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0 ) + , wellVariables_( wells_ ? (wells_arg->number_of_wells * numWellEq) : 0) + , F0_(wells_ ? (wells_arg->number_of_wells * numWellEq) : 0 ) { if( wells_ ) { @@ -32,9 +32,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, const VFPProperties* vfp_properties_arg, @@ -117,13 +117,9 @@ namespace Opm { - - - - template - template + template SimulatorReport - StandardWellsDense:: + StandardWellsDense:: assemble(Simulator& ebosSimulator, const int iterationIdx, const double dt, @@ -162,17 +158,16 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, bool only_wells) { - const int np = wells().number_of_phases; const int nw = wells().number_of_wells; + const int numComp = numComponents(); // clear all entries duneB_ = 0.0; @@ -191,45 +186,44 @@ namespace Opm { const int cell_idx = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - std::vector cq_s(np,0.0); - - std::vector mob(np, 0.0); + std::vector cq_s(numComp,0.0); + std::vector mob(numComp, 0.0); getMobility(ebosSimulator, perf, cell_idx, mob); - computeWellFlux(w, wells().WI[perf], intQuants.fluidState(), mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - for (int p1 = 0; p1 < np; ++p1) { + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. - const EvalWell cq_s_effective = cq_s[p1] * well_perforation_efficiency_factors_[perf]; + const EvalWell cq_s_effective = cq_s[componentIdx] * well_perforation_efficiency_factors_[perf]; if (!only_wells) { - // subtract sum of phase fluxes in the reservoir equation. + // subtract sum of component fluxes in the reservoir equation. // need to consider the efficiency factor - ebosResid[cell_idx][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.value(); + ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value(); } // subtract sum of phase fluxes in the well equations. - resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value(); + resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s[componentIdx].value(); // assemble the jacobians - for (int p2 = 0; p2 < np; ++p2) { + for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { if (!only_wells) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); - duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s_effective.derivative(p2+blocksize); // intput in transformed matrix - duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s_effective.derivative(p2); + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + duneB_[w][cell_idx][flowToEbosPvIdx(pvIdx)][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix + duneC_[w][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); } - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+blocksize); + invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s[componentIdx].derivative(pvIdx+numEq); } // add trivial equation for 2p cases (Only support water + oil) - if (np == 2) { + if (numComp == 2) { assert(!active_[ Gas ]); invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0; } // Store the perforation phase flux for later usage. - well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value(); + well_state.perfPhaseRates()[perf*numComp + componentIdx] = cq_s[componentIdx].value(); } // Store the perforation pressure for later usage. @@ -237,13 +231,13 @@ namespace Opm { } // add vol * dF/dt + Q to the well equations; - for (int p1 = 0; p1 < np; ++p1) { - EvalWell resWell_loc = (wellSurfaceVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; - resWell_loc += getQs(w, p1); - for (int p2 = 0; p2 < np; ++p2) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+blocksize); + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + EvalWell resWell_loc = (wellSurfaceVolumeFraction(w, componentIdx) - F0_[w + nw*componentIdx]) * volume / dt; + resWell_loc += getQs(w, componentIdx); + for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] += resWell_loc.derivative(pvIdx+numEq); } - resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value(); + resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] += resWell_loc.value(); } } @@ -252,15 +246,14 @@ namespace Opm { } - template - template + template void - StandardWellsDense:: + StandardWellsDense:: getMobility(const Simulator& ebosSimulator, const int perf, const int cell_idx, std::vector& mob) const { const int np = wells().number_of_phases; - assert (int(mob.size()) == np); + assert (int(mob.size()) == numComponents()); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); @@ -293,11 +286,11 @@ namespace Opm { - template - template + + template bool - StandardWellsDense:: - allow_cross_flow(const int w, Simulator& ebosSimulator) const + StandardWellsDense:: + allow_cross_flow(const int w, const Simulator& ebosSimulator) const { if (wells().allow_cf[w]) { return true; @@ -331,9 +324,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: localInvert(Mat& istlA) const { for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { @@ -348,9 +341,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: print(Mat& istlA) const { for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { @@ -364,9 +357,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: apply( BVector& r) const { if ( ! localWellsActive() ) { @@ -383,9 +376,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: apply(const BVector& x, BVector& Ax) { if ( ! localWellsActive() ) { @@ -406,9 +399,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) { if ( ! localWellsActive() ) { @@ -428,9 +421,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: recoverVariable(const BVector& x, BVector& xw) const { if ( ! localWellsActive() ) { @@ -445,11 +438,12 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { + assert(phaseIdx < 3); const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; return phaseToComp[ phaseIdx ]; } @@ -458,11 +452,12 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowToEbosPvIdx( const int flowPv ) const { + assert(flowPv < 3); const int flowToEbos[ 3 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, @@ -475,11 +470,12 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { + assert(phaseIdx < 3); const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; return flowToEbos[ phaseIdx ]; } @@ -488,11 +484,12 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: ebosCompToFlowPhaseIdx( const int compIdx ) const { + assert(compIdx < 3); const int compToPhase[ 3 ] = { Oil, Water, Gas }; return compToPhase[ compIdx ]; } @@ -501,9 +498,9 @@ namespace Opm { - template + template std::vector - StandardWellsDense:: + StandardWellsDense:: extractPerfData(const std::vector& in) const { const int nw = wells().number_of_wells; @@ -522,9 +519,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numPhases() const { return wells().number_of_phases; @@ -534,9 +531,9 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numCells() const { return pv_.size(); @@ -546,9 +543,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: resetWellControlFromState(const WellState& xw) const { const int nw = wells_->number_of_wells; @@ -562,9 +559,9 @@ namespace Opm { - template + template const Wells& - StandardWellsDense:: + StandardWellsDense:: wells() const { assert(wells_ != 0); @@ -575,9 +572,9 @@ namespace Opm { - template + template const Wells* - StandardWellsDense:: + StandardWellsDense:: wellsPointer() const { return wells_; @@ -587,9 +584,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: wellsActive() const { return wells_active_; @@ -599,9 +596,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: setWellsActive(const bool wells_active) { wells_active_ = wells_active; @@ -611,9 +608,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; @@ -623,27 +620,26 @@ namespace Opm { - template + template int - StandardWellsDense:: + StandardWellsDense:: numWellVars() const { if ( !localWellsActive() ) { return 0; } - // For each well, we have a bhp variable, and one flux per phase. const int nw = wells().number_of_wells; - return (numPhases() + 1) * nw; + return numWellEq * nw; } - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerforationDensities() const { return well_perforation_densities_; @@ -653,9 +649,9 @@ namespace Opm { - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; @@ -665,14 +661,14 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: extendEval(Eval in) const { EvalWell out = 0.0; out.setValue(in.value()); - for(int i = 0; i < blocksize;++i) { - out.setDerivative(i, in.derivative(flowToEbosPvIdx(i))); + for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { + out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx))); } return out; } @@ -681,18 +677,17 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: setWellVariables(const WellState& xw) { - const int np = wells().number_of_phases; const int nw = wells().number_of_wells; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { for (int w = 0; w < nw; ++w) { - wellVariables_[w + nw*phaseIdx] = 0.0; - wellVariables_[w + nw*phaseIdx].setValue(xw.wellSolutions()[w + nw* phaseIdx]); - wellVariables_[w + nw*phaseIdx].setDerivative(blocksize + phaseIdx, 1.0); + wellVariables_[w + nw*eqIdx] = 0.0; + wellVariables_[w + nw*eqIdx].setValue(xw.wellSolutions()[w + nw* eqIdx]); + wellVariables_[w + nw*eqIdx].setDerivative(numWellEq + eqIdx, 1.0); } } } @@ -701,9 +696,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: print(EvalWell in) const { std::cout << in.value() << std::endl; @@ -716,16 +711,15 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeAccumWells() { - const int np = wells().number_of_phases; const int nw = wells().number_of_wells; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { for (int w = 0; w < nw; ++w) { - F0_[w + nw * phaseIdx] = wellSurfaceVolumeFraction(w, phaseIdx).value(); + F0_[w + nw * eqIdx] = wellSurfaceVolumeFraction(w, eqIdx).value(); } } } @@ -734,28 +728,28 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellFlux(const int& w, const double& Tw, - const FluidState& fs, + const IntensiveQuantities& intQuants, const std::vector& mob_perfcells_dense, const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector& cq_s) const { const Opm::PhaseUsage& pu = phase_usage_; const int np = wells().number_of_phases; - std::vector cmix_s(np,0.0); - for (int phase = 0; phase < np; ++phase) { - //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - cmix_s[phase] = wellSurfaceVolumeFraction(w, phase); + const int numComp = numComponents(); + std::vector cmix_s(numComp,0.0); + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + cmix_s[componentIdx] = wellSurfaceVolumeFraction(w, componentIdx); } + auto& fs = intQuants.fluidState(); EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); EvalWell rs = extendEval(fs.Rs()); EvalWell rv = extendEval(fs.Rv()); - std::vector b_perfcells_dense(np, 0.0); + std::vector b_perfcells_dense(numComp, 0.0); for (int phase = 0; phase < np; ++phase) { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); @@ -772,26 +766,21 @@ namespace Opm { return; } - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - const EvalWell cq_p = - Tw * (mob_perfcells_dense[phase] * drawdown); - cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + // compute component volumetric rates at standard conditions + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown); + cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p; } if (active_[Oil] && active_[Gas]) { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; - const EvalWell cq_psOil = cq_ps[oilpos]; - const EvalWell cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs * cq_psOil; - cq_ps[oilpos] += rv * cq_psGas; + const EvalWell cq_sOil = cq_s[oilpos]; + const EvalWell cq_sGas = cq_s[gaspos]; + cq_s[gaspos] += rs * cq_sOil; + cq_s[oilpos] += rv * cq_sGas; } - // map to ADB - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase]; - } } else { //Do nothing if crossflow is not allowed if (!allow_cf && wells().type[w] == PRODUCER) { @@ -800,8 +789,8 @@ namespace Opm { // Using total mobilities EvalWell total_mob_dense = mob_perfcells_dense[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob_dense += mob_perfcells_dense[phase]; + for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) { + total_mob_dense += mob_perfcells_dense[componentIdx]; } // injection perforations total volume rates @@ -848,8 +837,8 @@ namespace Opm { // injecting connections total volumerates at standard conditions EvalWell cqt_is = cqt_i/volumeRatio; //std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl; - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cmix_s[phase] * cqt_is; // * b_perfcells_dense[phase]; + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; // * b_perfcells_dense[phase]; } } } @@ -858,10 +847,9 @@ namespace Opm { - template - template + template SimulatorReport - StandardWellsDense:: + StandardWellsDense:: solveWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state) @@ -921,9 +909,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: printIf(const int c, const double x, const double y, const double eps, const std::string type) const { if (std::abs(x-y) > eps) { @@ -935,9 +923,9 @@ namespace Opm { - template + template std::vector - StandardWellsDense:: + StandardWellsDense:: residual() const { if( ! wellsActive() ) @@ -945,14 +933,14 @@ namespace Opm { return std::vector(); } - const int np = numPhases(); const int nw = wells().number_of_wells; - std::vector res(np*nw); - for( int p=0; p res(numComp*nw); + for( int compIdx = 0; compIdx < numComp; ++compIdx) { + const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx); + for (int wellIdx = 0; wellIdx < nw; ++wellIdx) { + int idx = wellIdx + nw*compIdx; + res[idx] = resWell_[ wellIdx ][ ebosCompIdx ]; } } return res; @@ -962,10 +950,9 @@ namespace Opm { - template - template + template bool - StandardWellsDense:: + StandardWellsDense:: getWellConvergence(Simulator& ebosSimulator, const int iteration) const { @@ -973,11 +960,13 @@ namespace Opm { typedef std::vector< Scalar > Vector; const int np = numPhases(); + const int numComp = numComponents(); + const double tol_wells = param_.tolerance_wells_; const double maxResidualAllowed = param_.max_residual_allowed_; - std::vector< Scalar > B_avg( np, Scalar() ); - std::vector< Scalar > maxNormWell(np, Scalar() ); + std::vector< Scalar > B_avg( numComp, Scalar() ); + std::vector< Scalar > maxNormWell(numComp, Scalar() ); auto& grid = ebosSimulator.gridManager().grid(); const auto& gridView = grid.leafGridView(); @@ -993,10 +982,10 @@ 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 phaseIdx = 0; phaseIdx < np; ++phaseIdx ) { - auto& B = B_avg[ idx ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + auto& B = B_avg[ phaseIdx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); B += 1 / fs.invB(ebosPhaseIdx).value(); } @@ -1010,25 +999,25 @@ namespace Opm { } auto res = residual(); - const int nw = res.size() / np; + const int nw = res.size() / numComp; - for ( int idx = 0; idx < np; ++idx ) + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) { for ( int w = 0; w < nw; ++w ) { - maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(res[nw*idx + w])); + maxNormWell[compIdx] = std::max(maxNormWell[compIdx], std::abs(res[nw*compIdx + w])); } } grid.comm().max(maxNormWell.data(), maxNormWell.size()); - Vector well_flux_residual(np); + Vector well_flux_residual(numComp); bool converged_Well = true; // Finish computation - for ( int idx = 0; idx < np; ++idx ) + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) { - well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; - converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + well_flux_residual[compIdx] = B_avg[compIdx] * maxNormWell[compIdx]; + converged_Well = converged_Well && (well_flux_residual[compIdx] < tol_wells); } // if one of the residuals is NaN, throw exception, so that the solver can be restarted @@ -1061,8 +1050,8 @@ 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 phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - ss << std::setw(11) << well_flux_residual[phaseIdx]; + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + ss << std::setw(11) << well_flux_residual[compIdx]; } ss.precision(oprec); ss.flags(oflags); @@ -1075,10 +1064,9 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellConnectionPressures(const Simulator& ebosSimulator, const WellState& xw) { @@ -1099,10 +1087,9 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, const WellState& xw, std::vector& b_perf, @@ -1112,10 +1099,10 @@ namespace Opm { { const int nperf = wells().well_connpos[wells().number_of_wells]; const int nw = wells().number_of_wells; + const int numComp = numComponents(); const PhaseUsage& pu = phase_usage_; - const int np = phase_usage_.num_phases; - b_perf.resize(nperf*np); - surf_dens_perf.resize(nperf*np); + b_perf.resize(nperf*numComp); + surf_dens_perf.resize(nperf*numComp); //rs and rv are only used if both oil and gas is present if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { @@ -1136,16 +1123,16 @@ namespace Opm { const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); if (pu.phase_used[BlackoilPhases::Aqua]) { - b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * numComp] = FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); } if (pu.phase_used[BlackoilPhases::Vapour]) { - const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); if (oilrate > 0) { @@ -1168,11 +1155,11 @@ namespace Opm { } if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; if (pu.phase_used[BlackoilPhases::Vapour]) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; const double gasrate = std::abs(xw.wellRates()[gaspos_well]); if (gasrate > 0) { const double oilrate = std::abs(xw.wellRates()[oilpos_well]); @@ -1192,7 +1179,7 @@ namespace Opm { // Surface density. for (int p = 0; p < pu.num_phases; ++p) { - surf_dens_perf[np*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); } } } @@ -1202,9 +1189,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellState(const BVector& dwells, WellState& well_state) const { @@ -1484,9 +1471,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellControls(WellState& xw) const { // Even if there no wells active locally, we cannot @@ -1592,9 +1579,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateListEconLimited(const Schedule& schedule, const int current_step, const Wells* wells_struct, @@ -1700,9 +1687,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellConnectionDensitesPressures(const WellState& xw, const std::vector& b_perf, const std::vector& rsmax_perf, @@ -1714,7 +1701,7 @@ namespace Opm { // Compute densities well_perforation_densities_ = WellDensitySegmented::computeConnectionDensities( - wells(), xw, phase_usage_, + wells(), phase_usage_, xw.perfPhaseRates(), b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // Compute pressure deltas @@ -1727,10 +1714,9 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellPotentials(const Simulator& ebosSimulator, const WellState& well_state, std::vector& well_potentials) const @@ -1791,10 +1777,9 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: prepareTimeStep(const Simulator& ebos_simulator, WellState& well_state) { @@ -1877,10 +1862,9 @@ namespace Opm { - - template + template WellCollection* - StandardWellsDense:: + StandardWellsDense:: wellCollection() const { return well_collection_; @@ -1889,9 +1873,9 @@ namespace Opm { - template + template const std::vector& - StandardWellsDense:: + StandardWellsDense:: wellPerfEfficiencyFactors() const { return well_perforation_efficiency_factors_; @@ -1901,9 +1885,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: calculateEfficiencyFactors() { if ( !localWellsActive() ) { @@ -1929,9 +1913,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellVoidageRates(const WellState& well_state, std::vector& well_voidage_rates, std::vector& voidage_conversion_coeffs) const @@ -1989,9 +1973,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: applyVREPGroupControl(WellState& well_state) const { if ( wellCollection()->havingVREPGroups() ) { @@ -2017,9 +2001,9 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: getBhp(const int wellIdx) const { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == BHP) { @@ -2073,21 +2057,22 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - getQs(const int wellIdx, const int phaseIdx) const + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + getQs(const int wellIdx, const int compIdx) const { EvalWell qs = 0.0; const WellControls* wc = wells().ctrls[wellIdx]; const int np = wells().number_of_phases; + assert(compIdx < np); const int nw = wells().number_of_wells; const double target_rate = well_controls_get_current_target(wc); // TODO: the formulation for the injectors decides it only work with single phase // surface rate injection control. Improvement will be required. if (wells().type[wellIdx] == INJECTOR) { - const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + const double comp_frac = wells().comp_frac[np*wellIdx + compIdx]; if (comp_frac == 0.0) { return qs; } @@ -2101,7 +2086,7 @@ namespace Opm { // Producers if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { - return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); + return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx, compIdx); } if (well_controls_get_current_type(wc) == SURFACE_RATE) { @@ -2132,7 +2117,7 @@ namespace Opm { assert(phase_under_control >= 0); - if (phaseIdx == phase_under_control) { + if (compIdx == phase_under_control) { qs.setValue(target_rate); return qs; } @@ -2142,7 +2127,7 @@ namespace Opm { if (wellVolumeFractionScaled(wellIdx, phase_under_control) < eps) { return qs; } - return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx, phase_under_control)); + return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaled(wellIdx, phase_under_control)); } // when it is a combined two phase rate limit, such like LRAT @@ -2154,16 +2139,16 @@ namespace Opm { combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p); } } - return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / combined_volume_fraction); + return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / combined_volume_fraction); } // TODO: three phase surface rate control is not tested yet if (num_phases_under_rate_control == 3) { - return target_rate * wellSurfaceVolumeFraction(wellIdx, phaseIdx); + return target_rate * wellSurfaceVolumeFraction(wellIdx, compIdx); } } else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { // ReservoirRate - return target_rate * wellVolumeFractionScaled(wellIdx, phaseIdx); + return target_rate * wellVolumeFractionScaled(wellIdx, compIdx); } else { OPM_THROW(std::logic_error, "Unknown control type for well " << wells().name[wellIdx]); } @@ -2176,17 +2161,17 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - wellVolumeFraction(const int wellIdx, const int phaseIdx) const + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + wellVolumeFraction(const int wellIdx, const int compIdx) const { const int nw = wells().number_of_wells; - if (phaseIdx == Water) { + if (compIdx == Water) { return wellVariables_[WFrac * nw + wellIdx]; } - if (phaseIdx == Gas) { + if (compIdx == Gas) { return wellVariables_[GFrac * nw + wellIdx]; } @@ -2206,57 +2191,60 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + wellVolumeFractionScaled(const int wellIdx, const int compIdx) const { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { const double* distr = well_controls_get_current_distr(wc); - if (distr[phaseIdx] > 0.) { - return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx]; + assert(compIdx < 3); + if (distr[compIdx] > 0.) { + return wellVolumeFraction(wellIdx, compIdx) / distr[compIdx]; } else { // TODO: not sure why return EvalWell(0.) causing problem here // Probably due to the wrong Jacobians. - return wellVolumeFraction(wellIdx, phaseIdx); + return wellVolumeFraction(wellIdx, compIdx); } } + assert(compIdx < 3); std::vector g = {1,1,0.01}; - return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); + return (wellVolumeFraction(wellIdx, compIdx) / g[compIdx]); } - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - wellSurfaceVolumeFraction(const int well_index, const int phase) const + template + typename StandardWellsDense::EvalWell + StandardWellsDense:: + wellSurfaceVolumeFraction(const int well_index, const int compIdx) const { EvalWell sum_volume_fraction_scaled = 0.; - const int np = wells().number_of_phases; - for (int p = 0; p < np; ++p) { - sum_volume_fraction_scaled += wellVolumeFractionScaled(well_index, p); + const int numComp = numComponents(); + for (int idx = 0; idx < numComp; ++idx) { + sum_volume_fraction_scaled += wellVolumeFractionScaled(well_index, idx); } assert(sum_volume_fraction_scaled.value() != 0.); - return wellVolumeFractionScaled(well_index, phase) / sum_volume_fraction_scaled; + return wellVolumeFractionScaled(well_index, compIdx) / sum_volume_fraction_scaled; } - template + template bool - StandardWellsDense:: + StandardWellsDense:: checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const int well_number) const { + const Opm::PhaseUsage& pu = phase_usage_; const int np = well_state.numPhases(); @@ -2301,9 +2289,9 @@ namespace Opm { - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: + template + typename StandardWellsDense::RatioCheckTuple + StandardWellsDense:: checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const WellMapEntryType& map_entry) const @@ -2359,9 +2347,9 @@ namespace Opm { - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: + template + typename StandardWellsDense::RatioCheckTuple + StandardWellsDense:: checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, const WellState& well_state, const WellMapEntryType& map_entry) const @@ -2439,9 +2427,9 @@ namespace Opm { - template + template void - StandardWellsDense:: + StandardWellsDense:: updateWellStateWithTarget(const WellControls* wc, const int current, const int well_index, @@ -2656,9 +2644,9 @@ namespace Opm { - template + template bool - StandardWellsDense:: + StandardWellsDense:: wellHasTHPConstraints(const int well_index) const { const WellControls* well_control = wells().ctrls[well_index]; @@ -2675,10 +2663,9 @@ namespace Opm { - template - template + template void - StandardWellsDense:: + StandardWellsDense:: computeWellRatesWithBhp(const Simulator& ebosSimulator, const EvalWell& bhp, const int well_index, @@ -2695,7 +2682,7 @@ namespace Opm { std::vector cq_s(np, 0.0); std::vector mob(np, 0.0); getMobility(ebosSimulator, perf, cell_index, mob); - computeWellFlux(well_index, wells().WI[perf], intQuants.fluidState(), mob, bhp, + computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); for(int p = 0; p < np; ++p) { @@ -2709,9 +2696,9 @@ namespace Opm { - template + template double - StandardWellsDense:: + StandardWellsDense:: mostStrictBhpFromBhpLimits(const int well_index) const { double bhp; @@ -2765,10 +2752,9 @@ namespace Opm { - template - template + template std::vector - StandardWellsDense:: + StandardWellsDense:: computeWellPotentialWithTHP(const Simulator& ebosSimulator, const int well_index, const double initial_bhp, // bhp from BHP constraints diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 7f60d23af..aa4fef673 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -315,7 +315,7 @@ namespace Opm // Compute densities std::vector 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]; diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp index 576fac3a1..46f53f33a 100644 --- a/opm/autodiff/WellDensitySegmented.cpp +++ b/opm/autodiff/WellDensitySegmented.cpp @@ -28,12 +28,10 @@ - - std::vector Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells, - const WellStateFullyImplicitBlackoil& wstate, const PhaseUsage& phase_usage, + const std::vector& perfComponentRates, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& 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 q_out_perf(nperf*np); + std::vector 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 mix(np); - std::vector x(np); - std::vector surf_dens(np); + std::vector mix(numComponents,0.0); + std::vector x(numComponents); + std::vector surf_dens(numComponents); std::vector 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. diff --git a/opm/autodiff/WellDensitySegmented.hpp b/opm/autodiff/WellDensitySegmented.hpp index cf9a9d544..1091a823d 100644 --- a/opm/autodiff/WellDensitySegmented.hpp +++ b/opm/autodiff/WellDensitySegmented.hpp @@ -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 computeConnectionDensities(const Wells& wells, - const WellStateFullyImplicitBlackoil& wstate, const PhaseUsage& phase_usage, + const std::vector& perfComponentRates, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& 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) diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index b6271542d..de4b945f2 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -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 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]; diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp index 470aa9e48..520dc43ce 100644 --- a/tests/test_welldensitysegmented.cpp +++ b/tests/test_welldensitysegmented.cpp @@ -103,7 +103,7 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas) std::vector cd = WellDensitySegmented::computeConnectionDensities( - *wells, wellstate, pu, + *wells, pu, rates, b_perf, rsmax_perf, rvmax_perf, surf_dens); std::vector dp = WellDensitySegmented::computeConnectionPressureDelta(