diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 222edb60e..6cfb634cb 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -91,6 +91,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_solventprops_ad.cpp tests/test_multisegmentwells.cpp tests/test_multiphaseupwind.cpp + tests/test_wellmodel.cpp # tests/test_thresholdpressure.cpp tests/test_wellswitchlogger.cpp tests/test_timer.cpp @@ -102,6 +103,7 @@ list (APPEND TEST_DATA_FILES tests/VFPPROD2 tests/msw.data tests/TESTTIMER.DATA + tests/TESTWELLMODEL.DATA ) @@ -239,6 +241,10 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellHelpers.hpp opm/autodiff/StandardWells.hpp opm/autodiff/StandardWells_impl.hpp + opm/autodiff/WellInterface.hpp + opm/autodiff/WellInterface_impl.hpp + opm/autodiff/StandardWell.hpp + opm/autodiff/StandardWell_impl.hpp opm/autodiff/StandardWellsDense.hpp opm/autodiff/StandardWellsSolvent.hpp opm/autodiff/StandardWellsSolvent_impl.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 83e80cdc6..e2bbfdb78 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -154,7 +154,7 @@ namespace Opm { /// \param[in] terminal_output request output to cout/cerr BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, - const StandardWellsDense& well_model, + StandardWellsDense& well_model, RateConverterType& rate_converter, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output @@ -284,10 +284,9 @@ namespace Opm { const int nc = AutoDiffGrid::numCells(grid_); const int nw = numWells(); BVector x(nc); - BVector xw(nw); try { - solveJacobianSystem(x, xw); + solveJacobianSystem(x); report.linear_solve_time += perfTimer.stop(); report.total_linear_iterations += linearIterationsLastSolve(); } @@ -322,7 +321,11 @@ namespace Opm { // Apply the update, with considering model-dependent limitations and // chopping of the update. updateState(x,iteration); - wellModel().updateWellState(xw, well_state); + + if( nw > 0 ) + { + wellModel().recoverWellSolutionAndUpdateWellState(x, well_state); + } report.update_time += perfTimer.stop(); } @@ -488,7 +491,7 @@ namespace Opm { /// Solve the Jacobian system Jx = r where J is the Jacobian and /// r is the residual. - void solveJacobianSystem(BVector& x, BVector& xw) const + void solveJacobianSystem(BVector& x) const { const auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); @@ -510,13 +513,6 @@ namespace Opm { Operator opA(ebosJac, well_model_); istlSolver().solve( opA, x, ebosResid ); } - - if( xw.size() > 0 ) - { - // recover wells. - xw = 0.0; - wellModel().recoverVariable(x, xw); - } } //===================================================================== @@ -766,8 +762,7 @@ namespace Opm { const double pvSumLocal, std::vector< Scalar >& R_sum, std::vector< Scalar >& maxCoeff, - std::vector< Scalar >& B_avg, - std::vector< Scalar >& maxNormWell ) + std::vector< Scalar >& B_avg) { // Compute total pore volume (use only owned entries) double pvSum = pvSumLocal; @@ -779,13 +774,12 @@ namespace Opm { std::vector< Scalar > maxBuffer; const int numComp = B_avg.size(); sumBuffer.reserve( 2*numComp + 1 ); // +1 for pvSum - maxBuffer.reserve( 2*numComp ); + maxBuffer.reserve( numComp ); for( int compIdx = 0; compIdx < numComp; ++compIdx ) { 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 @@ -801,11 +795,14 @@ namespace Opm { for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx ) { B_avg[ compIdx ] = sumBuffer[ buffIdx ]; - maxCoeff[ compIdx ] = maxBuffer[ buffIdx ]; ++buffIdx; R_sum[ compIdx ] = sumBuffer[ buffIdx ]; - maxNormWell[ compIdx ] = maxBuffer[ buffIdx ]; + } + + for( int compIdx = 0; compIdx < numComp; ++compIdx ) + { + maxCoeff[ compIdx ] = maxBuffer[ compIdx ]; } // restore global pore volume @@ -828,7 +825,6 @@ namespace Opm { const double dt = timer.currentStepLength(); const double tol_mb = param_.tolerance_mb_; const double tol_cnv = param_.tolerance_cnv_; - const double tol_wells = param_.tolerance_wells_; const int np = numPhases(); const int numComp = numComponents(); @@ -836,7 +832,6 @@ namespace Opm { Vector R_sum(numComp, 0.0 ); Vector B_avg(numComp, 0.0 ); Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); - Vector maxNormWell(numComp, 0.0 ); const auto& ebosModel = ebosSimulator_.model(); const auto& ebosProblem = ebosSimulator_.problem(); @@ -896,24 +891,16 @@ namespace Opm { B_avg[ i ] /= Scalar( global_nc_ ); } - // compute maximum of local well residuals - const Vector& wellResidual = wellModel().residual(); - const int nw = wellResidual.size() / numComp; - assert(nw * numComp == int(wellResidual.size())); - for( int compIdx = 0; compIdx < numComp; ++compIdx ) - { - for ( int w = 0; w < nw; ++w ) { - maxNormWell[compIdx] = std::max(maxNormWell[compIdx], std::abs(wellResidual[nw*compIdx + w])); - } - } + // TODO: we remove the maxNormWell for now because the convergence of wells are on a individual well basis. + // Anyway, we need to provide some infromation to help debug the well iteration process. + // compute global sum and max of quantities const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal, - R_sum, maxCoeff, B_avg, maxNormWell ); + R_sum, maxCoeff, B_avg); Vector CNV(numComp); Vector mass_balance_residual(numComp); - Vector well_flux_residual(numComp); bool converged_MB = true; bool converged_CNV = true; @@ -927,8 +914,7 @@ namespace Opm { converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv); // Well flux convergence is only for fluid phases, not other materials // in our current implementation. - well_flux_residual[compIdx] = B_avg[compIdx] * maxNormWell[compIdx]; - converged_Well = converged_Well && (well_flux_residual[compIdx] < tol_wells); + converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg); residual_norms.push_back(CNV[compIdx]); } @@ -965,9 +951,6 @@ namespace Opm { for (int compIdx = 0; compIdx < numComp; ++compIdx) { msg += " CNV(" + key[ compIdx ] + ") "; } - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - msg += " W-FLUX(" + key[ compIdx ] + ")"; - } OpmLog::note(msg); } std::ostringstream ss; @@ -980,9 +963,6 @@ namespace Opm { for (int compIdx = 0; compIdx < numComp; ++compIdx) { ss << std::setw(11) << CNV[compIdx]; } - 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()); @@ -992,13 +972,11 @@ namespace Opm { const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); if (std::isnan(mass_balance_residual[phaseIdx]) - || std::isnan(CNV[phaseIdx]) - || (phaseIdx < numPhases() && std::isnan(well_flux_residual[phaseIdx]))) { + || std::isnan(CNV[phaseIdx])) { OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); } if (mass_balance_residual[phaseIdx] > maxResidualAllowed() - || CNV[phaseIdx] > maxResidualAllowed() - || (phaseIdx < numPhases() && well_flux_residual[phaseIdx] > maxResidualAllowed())) { + || CNV[phaseIdx] > maxResidualAllowed()) { OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); } } @@ -1523,7 +1501,7 @@ namespace Opm { SimulatorReport failureReport_; // Well Model - StandardWellsDense well_model_; + StandardWellsDense& well_model_; /// \brief Whether we print something to std::cout bool terminal_output_; @@ -1545,13 +1523,7 @@ namespace Opm { const StandardWellsDense& wellModel() const { return well_model_; } - /// return the Well struct in the StandardWells - const Wells& wells() const { return well_model_.wells(); } - - /// return true if wells are available in the reservoir - bool wellsActive() const { return well_model_.wellsActive(); } - - int numWells() const { return wellsActive() ? wells().number_of_wells : 0; } + int numWells() const { return well_model_.numWells(); } /// return true if wells are available on this process bool localWellsActive() const { return well_model_.localWellsActive(); } diff --git a/opm/autodiff/Compat.cpp b/opm/autodiff/Compat.cpp index e346f4ad4..289255df4 100644 --- a/opm/autodiff/Compat.cpp +++ b/opm/autodiff/Compat.cpp @@ -294,9 +294,6 @@ void wellsToState( const data::Wells& wells, { // Set base class variables. wellsToState(wells, phases, static_cast(state)); - - // Set wellSolution() variable. - state.setWellSolutions(phases); } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index cc7b4bc6a..5085427fa 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -141,7 +141,6 @@ public: ExtraData extra; failureReport_ = SimulatorReport(); - extractLegacyPoreVolume_(); extractLegacyDepth_(); // communicate the initial solution to ebos @@ -480,7 +479,6 @@ protected: activePhases, gravity, legacyDepth_, - legacyPoreVolume_, globalNumCells, grid()); auto model = std::unique_ptr(new Model(ebosSimulator_, @@ -867,22 +865,6 @@ protected: } } - void extractLegacyPoreVolume_() - { - const auto& grid = ebosSimulator_.gridManager().grid(); - const unsigned numCells = grid.size(/*codim=*/0); - const auto& ebosProblem = ebosSimulator_.problem(); - const auto& ebosModel = ebosSimulator_.model(); - - legacyPoreVolume_.resize(numCells); - for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) { - // todo (?): respect rock compressibility - legacyPoreVolume_[cellIdx] = - ebosModel.dofTotalVolume(cellIdx) - *ebosProblem.porosity(cellIdx); - } - } - void extractLegacyDepth_() { const auto& grid = ebosSimulator_.gridManager().grid(); @@ -1009,7 +991,6 @@ protected: Simulator& ebosSimulator_; std::vector legacyCellPvtRegionIdx_; - std::vector legacyPoreVolume_; std::vector legacyDepth_; typedef typename Solver::SolverParameters SolverParameters; diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp new file mode 100644 index 000000000..036ba51c1 --- /dev/null +++ b/opm/autodiff/StandardWell.hpp @@ -0,0 +1,305 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + Copyright 2016 - 2017 IRIS AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_STANDARDWELL_HEADER_INCLUDED +#define OPM_STANDARDWELL_HEADER_INCLUDED + + +#include + +namespace Opm +{ + + template + class StandardWell: public WellInterface + { + + public: + typedef WellInterface Base; + // TODO: some functions working with AD variables handles only with values (double) without + // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value. + // And also, it can also be beneficial to make these functions hanle different types of AD variables. + using typename Base::Simulator; + using typename Base::WellState; + using typename Base::IntensiveQuantities; + using typename Base::FluidSystem; + using typename Base::MaterialLaw; + using typename Base::ModelParameters; + using typename Base::BlackoilIndices; + using typename Base::PolymerModule; + + // the positions of the primary variables for StandardWell + // there are three primary variables, the second and the third ones are F_w and F_g + // the first one can be total rate (G_t) or bhp, based on the control + enum WellVariablePositions { + XvarWell = 0, + WFrac = 1, + GFrac = 2, + SFrac = 3 + }; + + using typename Base::Scalar; + using typename Base::ConvergenceReport; + + using Base::numEq; + + using Base::has_solvent; + using Base::has_polymer; + using Base::name; + + // TODO: with flow_ebos，for a 2P deck, // TODO: for the 2p deck, numEq will be 3, a dummy phase is already added from the reservoir side. + // it will cause problem here without processing the dummy phase. + static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // number of wellEq is only numEq - 1 for polymer + using typename Base::Mat; + using typename Base::BVector; + using typename Base::Eval; + + // sparsity pattern for the matrices + //[A C^T [x = [ res + // B D ] x_well] res_well] + + // the vector type for the res_well and x_well + typedef Dune::FieldVector VectorBlockWellType; + typedef Dune::BlockVector BVectorWell; + + // the matrix type for the diagonal matrix D + typedef Dune::FieldMatrix DiagMatrixBlockWellType; + typedef Dune::BCRSMatrix DiagMatWell; + + // the matrix type for the non-diagonal matrix B and C^T + typedef Dune::FieldMatrix OffDiagMatrixBlockWellType; + typedef Dune::BCRSMatrix OffDiagMatWell; + + typedef DenseAd::Evaluation EvalWell; + + // TODO: should these go to WellInterface? + static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; + + + StandardWell(const Well* well, const int time_step, const Wells* wells); + + virtual void init(const PhaseUsage* phase_usage_arg, + const std::vector* active_arg, + const std::vector& depth_arg, + const double gravity_arg, + const int num_cells); + + + virtual void initPrimaryVariablesEvaluation() const; + + virtual void assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells); + + /// updating the well state based the control mode specified with current + // TODO: later will check wheter we need current + virtual void updateWellStateWithTarget(const int current, + WellState& xw) const; + + // TODO: this should go to the WellInterface, while updateWellStateWithTarget + // will need touch different types of well_state, we will see. + virtual void updateWellControl(WellState& xw, + wellhelpers::WellSwitchingLogger& logger) const; + + /// check whether the well equations get converged for this well + virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg, + const ModelParameters& param) const; + + /// computing the accumulation term for later use in well mass equations + virtual void computeAccumWell(); + + virtual void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw); + + /// Ax = Ax - C D^-1 B x + virtual void apply(const BVector& x, BVector& Ax) const; + /// r = r - C D^-1 Rw + virtual void apply(BVector& r) const; + + /// using the solution x to recover the solution xw for wells and applying + /// xw to update Well State + virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param, + WellState& well_state) const; + + /// computing the well potentials for group control + virtual void computeWellPotentials(const Simulator& ebosSimulator, + const WellState& well_state, + std::vector& well_potentials) const; + + virtual void updatePrimaryVariables(const WellState& well_state) const; + + protected: + + // protected functions from the Base class + using Base::getAllowCrossFlow; + using Base::phaseUsage; + using Base::active; + using Base::flowToEbosPvIdx; + using Base::flowPhaseToEbosPhaseIdx; + using Base::flowPhaseToEbosCompIdx; + using Base::numComponents; + using Base::wsolvent; + using Base::wpolymer; + using Base::wellHasTHPConstraints; + using Base::mostStrictBhpFromBhpLimits; + + // protected member variables from the Base class + using Base::vfp_properties_; + using Base::gravity_; + using Base::well_efficiency_factor_; + using Base::phase_usage_; + using Base::first_perf_; + using Base::ref_depth_; + using Base::perf_depth_; + using Base::well_cells_; + using Base::number_of_perforations_; + using Base::number_of_phases_; + using Base::saturation_table_number_; + using Base::comp_frac_; + using Base::well_index_; + using Base::index_of_well_; + using Base::well_controls_; + using Base::well_type_; + + using Base::perf_rep_radius_; + using Base::perf_length_; + using Base::bore_diameters_; + + // densities of the fluid in each perforation + std::vector perf_densities_; + // pressure drop between different perforations + std::vector perf_pressure_diffs_; + + // residuals of the well equations + BVectorWell resWell_; + + // two off-diagonal matrices + OffDiagMatWell duneB_; + OffDiagMatWell duneC_; + // diagonal matrix for the well + DiagMatWell invDuneD_; + + // several vector used in the matrix calculation + mutable BVectorWell Bx_; + mutable BVectorWell invDrw_; + + // the values for the primary varibles + // based on different solutioin strategies, the wells can have different primary variables + mutable std::vector primary_variables_; + + // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation + mutable std::vector primary_variables_evaluation_; + + // the saturations in the well bore under surface conditions at the beginning of the time step + std::vector F0_; + + // TODO: this function should be moved to the base class. + // while it faces chanllenges for MSWell later, since the calculation of bhp + // based on THP is never implemented for MSWell yet. + EvalWell getBhp() const; + + // TODO: it is also possible to be moved to the base class. + EvalWell getQs(const int comp_idx) const; + + EvalWell wellVolumeFractionScaled(const int phase) const; + + EvalWell wellVolumeFraction(const int phase) const; + + EvalWell wellSurfaceVolumeFraction(const int phase) const; + + EvalWell extendEval(const Eval& in) const; + + bool crossFlowAllowed(const Simulator& ebosSimulator) const; + + // xw = inv(D)*(rw - C*x) + void recoverSolutionWell(const BVector& x, BVectorWell& xw) const; + + // updating the well_state based on well solution dwells + void updateWellState(const BVectorWell& dwells, + const BlackoilModelParameters& param, + WellState& well_state) const; + + // calculate the properties for the well connections + // to calulate the pressure difference between well connections. + 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; + + // TODO: not total sure whether it is a good idea to put this function here + // the major reason to put here is to avoid the usage of Wells struct + void computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf); + + void computeConnectionPressureDelta(); + + void computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf); + + virtual void solveEqAndUpdateWellState(const ModelParameters& param, + WellState& well_state); + + // TODO: to check whether all the paramters are required + void computePerfRate(const IntensiveQuantities& intQuants, + const std::vector& mob_perfcells_dense, + const double Tw, const EvalWell& bhp, const double& cdp, + const bool& allow_cf, std::vector& cq_s) const; + + // TODO: maybe we should provide a light version of computePerfRate, which does not include the + // calculation of the derivatives + void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const EvalWell& bhp, + std::vector& well_flux) const; + + std::vector computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const double initial_bhp, // bhp from BHP constraints + const std::vector& initial_potential) const; + + template + ValueType calculateBhpFromThp(const std::vector& rates, const int control_index) const; + + double calculateThpFromBhp(const std::vector& rates, const int control_index, const double bhp) const; + + // get the mobility for specific perforation + void getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const; + }; + +} + +#include "StandardWell_impl.hpp" + +#endif // OPM_STANDARDWELL_HEADER_INCLUDED diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp new file mode 100644 index 000000000..70f8782fa --- /dev/null +++ b/opm/autodiff/StandardWell_impl.hpp @@ -0,0 +1,2060 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + Copyright 2016 - 2017 IRIS AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +namespace Opm +{ + template + StandardWell:: + StandardWell(const Well* well, const int time_step, const Wells* wells) + : Base(well, time_step, wells) + , perf_densities_(number_of_perforations_) + , perf_pressure_diffs_(number_of_perforations_) + , primary_variables_(numWellEq, 0.0) + , primary_variables_evaluation_(numWellEq) // the number of the primary variables + , F0_(numWellEq) + { + duneB_.setBuildMode( OffDiagMatWell::row_wise ); + duneC_.setBuildMode( OffDiagMatWell::row_wise ); + invDuneD_.setBuildMode( DiagMatWell::row_wise ); + } + + + + + + template + void + StandardWell:: + init(const PhaseUsage* phase_usage_arg, + const std::vector* active_arg, + const std::vector& depth_arg, + const double gravity_arg, + const int num_cells) + { + Base::init(phase_usage_arg, active_arg, + depth_arg, gravity_arg, num_cells); + + perf_depth_.resize(number_of_perforations_, 0.); + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + perf_depth_[perf] = depth_arg[cell_idx]; + } + + // setup sparsity pattern for the matrices + //[A C^T [x = [ res + // B D] x_well] res_well] + // set the size of the matrices + invDuneD_.setSize(1, 1, 1); + duneB_.setSize(1, num_cells, number_of_perforations_); + duneC_.setSize(1, num_cells, number_of_perforations_); + + for (auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + row.insert(row.index()); + } + + for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + for (int perf = 0 ; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + row.insert(cell_idx); + } + } + + // make the C^T matrix + for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) { + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + row.insert(cell_idx); + } + } + + resWell_.resize(1); + + // resize temporary class variables + Bx_.resize( duneB_.N() ); + invDrw_.resize( invDuneD_.N() ); + } + + + + + + template + void StandardWell:: + initPrimaryVariablesEvaluation() const + { + // TODO: using numComp here is only to make the 2p + dummy phase work + // TODO: in theory, we should use numWellEq here. + // for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { + for (int eqIdx = 0; eqIdx < numComponents(); ++eqIdx) { + assert( eqIdx < primary_variables_.size() ); + + primary_variables_evaluation_[eqIdx] = 0.0; + primary_variables_evaluation_[eqIdx].setValue(primary_variables_[eqIdx]); + primary_variables_evaluation_[eqIdx].setDerivative(numEq + eqIdx, 1.0); + } + } + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + getBhp() const + { + const WellControls* wc = well_controls_; + if (well_controls_get_current_type(wc) == BHP) { + EvalWell bhp = 0.0; + const double target_rate = well_controls_get_current_target(wc); + bhp.setValue(target_rate); + return bhp; + } else if (well_controls_get_current_type(wc) == THP) { + const int control = well_controls_get_current(wc); + + const Opm::PhaseUsage& pu = phaseUsage(); + std::vector rates(3, 0.0); + if (active()[ Water ]) { + rates[ Water ]= getQs(pu.phase_pos[ Water]); + } + if (active()[ Oil ]) { + rates[ Oil ] = getQs(pu.phase_pos[ Oil ]); + } + if (active()[ Gas ]) { + rates[ Gas ] = getQs(pu.phase_pos[ Gas ]); + } + return calculateBhpFromThp(rates, control); + } + + return primary_variables_evaluation_[XvarWell]; + } + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + getQs(const int comp_idx) const // TODO: phase or component? + { + EvalWell qs = 0.0; + + const WellControls* wc = well_controls_; + const int np = number_of_phases_; + const double target_rate = well_controls_get_current_target(wc); + + assert(comp_idx < numComponents()); + const auto pu = phaseUsage(); + + // TODO: the formulation for the injectors decides it only work with single phase + // surface rate injection control. Improvement will be required. + if (well_type_ == INJECTOR) { + if (has_solvent) { + // TODO: investigate whether the use of the comp_frac is justified. + // The usage of the comp_frac is not correct, which should be changed later. + double comp_frac = 0.0; + if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent + comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent(); + } else if (comp_idx == pu.phase_pos[ Gas ]) { + comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent()); + } else { + comp_frac = comp_frac_[comp_idx]; + } + if (comp_frac == 0.0) { + return qs; //zero + } + + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + return comp_frac * primary_variables_evaluation_[XvarWell]; + } + + qs.setValue(comp_frac * target_rate); + return qs; + } + + const double comp_frac = comp_frac_[comp_idx]; + if (comp_frac == 0.0) { + return qs; + } + + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + return primary_variables_evaluation_[XvarWell]; + } + qs.setValue(target_rate); + return qs; + } + + // Producers + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { + return primary_variables_evaluation_[XvarWell] * wellVolumeFractionScaled(comp_idx); + } + + if (well_controls_get_current_type(wc) == SURFACE_RATE) { + // checking how many phases are included in the rate control + // to decide wheter it is a single phase rate control or not + const double* distr = well_controls_get_current_distr(wc); + int num_phases_under_rate_control = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + num_phases_under_rate_control += 1; + } + } + + // there should be at least one phase involved + assert(num_phases_under_rate_control > 0); + + // when it is a single phase rate limit + if (num_phases_under_rate_control == 1) { + + // looking for the phase under control + int phase_under_control = -1; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + phase_under_control = phase; + break; + } + } + + assert(phase_under_control >= 0); + + EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(phase_under_control); + if (has_solvent && phase_under_control == Gas) { + // for GRAT controlled wells solvent is included in the target + wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx); + } + + if (comp_idx == phase_under_control) { + if (has_solvent && phase_under_control == Gas) { + qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); + return qs; + } + qs.setValue(target_rate); + return qs; + } + + // TODO: not sure why the single phase under control will have near zero fraction + const double eps = 1e-6; + if (wellVolumeFractionScaledPhaseUnderControl < eps) { + return qs; + } + return (target_rate * wellVolumeFractionScaled(comp_idx) / wellVolumeFractionScaledPhaseUnderControl); + } + + // when it is a combined two phase rate limit, such like LRAT + // we neec to calculate the rate for the certain phase + if (num_phases_under_rate_control == 2) { + EvalWell combined_volume_fraction = 0.; + for (int p = 0; p < np; ++p) { + if (distr[p] == 1.0) { + combined_volume_fraction += wellVolumeFractionScaled(p); + } + } + return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction); + } + + // TODO: three phase surface rate control is not tested yet + if (num_phases_under_rate_control == 3) { + return target_rate * wellSurfaceVolumeFraction(comp_idx); + } + } else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + // ReservoirRate + return target_rate * wellVolumeFractionScaled(comp_idx); + } else { + OPM_THROW(std::logic_error, "Unknown control type for well " << name()); + } + + // avoid warning of condition reaches end of non-void function + return qs; + } + + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + wellVolumeFractionScaled(const int compIdx) const + { + // TODO: we should be able to set the g for the well based on the control type + // instead of using explicit code for g all the times + const WellControls* wc = well_controls_; + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + + if (has_solvent && compIdx == contiSolventEqIdx) { + return wellVolumeFraction(compIdx); + } + const double* distr = well_controls_get_current_distr(wc); + assert(compIdx < 3); + if (distr[compIdx] > 0.) { + return wellVolumeFraction(compIdx) / distr[compIdx]; + } else { + // TODO: not sure why return EvalWell(0.) causing problem here + // Probably due to the wrong Jacobians. + return wellVolumeFraction(compIdx); + } + } + + std::vector g = {1, 1, 0.01, 0.01}; + return (wellVolumeFraction(compIdx) / g[compIdx]); + } + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + wellVolumeFraction(const int compIdx) const + { + if (compIdx == Water) { + return primary_variables_evaluation_[WFrac]; + } + + if (compIdx == Gas) { + return primary_variables_evaluation_[GFrac]; + } + + if (has_solvent && compIdx == contiSolventEqIdx) { + return primary_variables_evaluation_[SFrac]; + } + + // Oil fraction + EvalWell well_fraction = 1.0; + if (active()[Water]) { + well_fraction -= primary_variables_evaluation_[WFrac]; + } + + if (active()[Gas]) { + well_fraction -= primary_variables_evaluation_[GFrac]; + } + if (has_solvent) { + well_fraction -= primary_variables_evaluation_[SFrac]; + } + return well_fraction; + } + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + wellSurfaceVolumeFraction(const int compIdx) const + { + EvalWell sum_volume_fraction_scaled = 0.; + const int numComp = numComponents(); + for (int idx = 0; idx < numComp; ++idx) { + sum_volume_fraction_scaled += wellVolumeFractionScaled(idx); + } + + assert(sum_volume_fraction_scaled.value() != 0.); + + return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled; + } + + + + + + template + typename StandardWell::EvalWell + StandardWell:: + extendEval(const Eval& in) const + { + EvalWell out = 0.0; + out.setValue(in.value()); + for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { + out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx))); + } + return out; + } + + + + + + template + void + StandardWell:: + computePerfRate(const IntensiveQuantities& intQuants, + const std::vector& mob_perfcells_dense, + const double Tw, const EvalWell& bhp, const double& cdp, + const bool& allow_cf, std::vector& cq_s) const + { + const Opm::PhaseUsage& pu = phaseUsage(); + const int np = number_of_phases_; + const int numComp = numComponents(); + std::vector cmix_s(numComp,0.0); + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + cmix_s[componentIdx] = wellSurfaceVolumeFraction(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(numComp, 0.0); + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + } + if (has_solvent) { + b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); + } + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + cdp; + EvalWell drawdown = pressure - well_pressure; + + // producing perforations + if ( drawdown.value() > 0 ) { + //Do nothing if crossflow is not allowed + if (!allow_cf && well_type_ == INJECTOR) { + return; + } + + // 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_sOil = cq_s[oilpos]; + const EvalWell cq_sGas = cq_s[gaspos]; + cq_s[gaspos] += rs * cq_sOil; + cq_s[oilpos] += rv * cq_sGas; + } + + } else { + //Do nothing if crossflow is not allowed + if (!allow_cf && well_type_ == PRODUCER) { + return; + } + + // Using total mobilities + EvalWell total_mob_dense = mob_perfcells_dense[0]; + for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) { + total_mob_dense += mob_perfcells_dense[componentIdx]; + } + + // injection perforations total volume rates + const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio = 0.0; + if (active()[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if (has_solvent) { + volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx]; + } + + if (active()[Oil] && active()[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + + // Incorporate RS/RV factors if both oil and gas active + const EvalWell d = 1.0 - rv * rs; + + if (d.value() == 0.0) { + OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << name() << " during flux calcuation" + << " with rs " << rs << " and rv " << rv); + } + + const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < + void + StandardWell:: + assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) + { + // TODO: accessing well_state information is the only place to use nw at the moment + const int nw = well_state.bhp().size(); + const int numComp = numComponents(); + const int np = number_of_phases_; + + // clear all entries + duneB_ = 0.0; + duneC_ = 0.0; + invDuneD_ = 0.0; + resWell_ = 0.0; + + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + // TODO: it probably can be static member for StandardWell + const double volume = 0.002831684659200; // 0.1 cu ft; + + const bool allow_cf = crossFlowAllowed(ebosSimulator); + + const EvalWell& bhp = getBhp(); + + for (int perf = 0; perf < number_of_perforations_; ++perf) { + + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector cq_s(numComp,0.0); + std::vector mob(numComp, 0.0); + getMobility(ebosSimulator, perf, mob); + computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + + 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[componentIdx] * well_efficiency_factor_; + + if (!only_wells) { + // subtract sum of component fluxes in the reservoir equation. + // need to consider the efficiency factor + ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value(); + } + + // subtract sum of phase fluxes in the well equations. + resWell_[0][componentIdx] -= cq_s[componentIdx].value(); + + // assemble the jacobians + for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + duneC_[0][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix + } + invDuneD_[0][0][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); + } + + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + duneB_[0][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + } + } + + // add a trivial equation for the dummy phase for 2p cases (Only support water + oil) + if ( numComp < numWellEq ) { + assert(!active()[ Gas ]); + invDuneD_[0][0][Gas][Gas] = 1.0; + } + + // Store the perforation phase flux for later usage. + if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) + well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); + } else { + well_state.perfPhaseRates()[(first_perf_ + perf) * np + componentIdx] = cq_s[componentIdx].value(); + } + } + + if (has_polymer) { + EvalWell cq_s_poly = cq_s[Water]; + if (well_type_ == INJECTOR) { + cq_s_poly *= wpolymer(); + } else { + cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); + } + if (!only_wells) { + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { + ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx); + } + ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); + } + } + + // Store the perforation pressure for later usage. + well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; + } + + // add vol * dF/dt + Q to the well equations; + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; + resWell_loc += getQs(componentIdx); + for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { + invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq); + } + resWell_[0][componentIdx] += resWell_loc.value(); + } + + // do the local inversion of D. + invDuneD_[0][0].invert(); + } + + + + + + template + bool + StandardWell:: + crossFlowAllowed(const Simulator& ebosSimulator) const + { + if (getAllowCrossFlow()) { + return true; + } + + // TODO: investigate the justification of the following situation + + // check for special case where all perforations have cross flow + // then the wells must allow for cross flow + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell bhp = getBhp(); + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + perf_pressure_diffs_[perf]; + EvalWell drawdown = pressure - well_pressure; + + if (drawdown.value() < 0 && well_type_ == INJECTOR) { + return false; + } + + if (drawdown.value() > 0 && well_type_ == PRODUCER) { + return false; + } + } + return true; + } + + + + + + template + void + StandardWell:: + getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const + { + const int np = number_of_phases_; + const int cell_idx = well_cells_[perf]; + assert (int(mob.size()) == numComponents()); + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + + // either use mobility of the perforation cell or calcualte its own + // based on passing the saturation table index + const int satid = saturation_table_number_[perf] - 1; + const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); + if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell + + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + } + if (has_solvent) { + mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); + } + } else { + + const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); + Eval relativePerms[3] = { 0.0, 0.0, 0.0 }; + MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); + + // reset the satnumvalue back to original + materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); + + // compute the mobility + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); + } + + // this may not work if viscosity and relperms has been modified? + if (has_solvent) { + OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent"); + } + } + + // modify the water mobility if polymer is present + if (has_polymer) { + // assume fully mixture for wells. + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + + if (well_type_ == INJECTOR) { + const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex()); + mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) ); + } + + if (PolymerModule::hasPlyshlog()) { + // compute the well water velocity with out shear effects. + const int numComp = numComponents(); + const bool allow_cf = crossFlowAllowed(ebosSimulator); + const EvalWell& bhp = getBhp(); + std::vector cq_s(numComp,0.0); + computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + // TODO: make area a member + double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf]; + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + const auto& scaledDrainageInfo = + materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx); + const Scalar& Swcr = scaledDrainageInfo.Swcr; + const EvalWell poro = extendEval(intQuants.porosity()); + const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); + // guard against zero porosity and no water + const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); + EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); + + if (PolymerModule::hasShrate()) { + // TODO Use the same conversion as for the reservoar equations. + // Need the "permeability" of the well? + // For now use the same formula as in legacy. + waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / bore_diameters_[perf]; + } + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, + intQuants.pvtRegionIndex(), + waterVelocity); + + // modify the mobility with the shear factor and recompute the well fluxes. + mob[ Water ] /= shearFactor; + } + } + } + + + + + + template + void + StandardWell:: + updateWellState(const BVectorWell& dwells, + const BlackoilModelParameters& param, + WellState& well_state) const + { + const int np = number_of_phases_; + const int nw = well_state.bhp().size(); + const double dBHPLimit = param.dbhp_max_rel_; + const double dFLimit = param.dwell_fraction_max_; + + const std::vector xvar_well_old = primary_variables_; + + // update the second and third well variable (The flux fractions) + std::vector F(np,0.0); + if (active()[ Water ]) { + const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit); + primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited; + } + + if (active()[ Gas ]) { + const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit); + primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited; + } + + if (has_solvent) { + const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]),dFLimit); + primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited; + } + + assert(active()[ Oil ]); + F[Oil] = 1.0; + + if (active()[ Water ]) { + F[Water] = primary_variables_[WFrac]; + F[Oil] -= F[Water]; + } + + if (active()[ Gas ]) { + F[Gas] = primary_variables_[GFrac]; + F[Oil] -= F[Gas]; + } + + double F_solvent = 0.0; + if (has_solvent) { + F_solvent = primary_variables_[SFrac]; + F[Oil] -= F_solvent; + } + + if (active()[ Water ]) { + if (F[Water] < 0.0) { + if (active()[ Gas ]) { + F[Gas] /= (1.0 - F[Water]); + } + if (has_solvent) { + F_solvent /= (1.0 - F[Water]); + } + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + } + + if (active()[ Gas ]) { + if (F[Gas] < 0.0) { + if (active()[ Water ]) { + F[Water] /= (1.0 - F[Gas]); + } + if (has_solvent) { + F_solvent /= (1.0 - F[Gas]); + } + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } + } + + if (F[Oil] < 0.0) { + if (active()[ Water ]) { + F[Water] /= (1.0 - F[Oil]); + } + if (active()[ Gas ]) { + F[Gas] /= (1.0 - F[Oil]); + } + if (has_solvent) { + F_solvent /= (1.0 - F[Oil]); + } + F[Oil] = 0.0; + } + + if (active()[ Water ]) { + primary_variables_[WFrac] = F[Water]; + } + if (active()[ Gas ]) { + primary_variables_[GFrac] = F[Gas]; + } + if(has_solvent) { + primary_variables_[SFrac] = F_solvent; + } + + // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. + // More testing is needed to make sure this is correct for well groups and THP. + if (has_solvent){ + F[Gas] += F_solvent; + } + + // The interpretation of the first well variable depends on the well control + const WellControls* wc = well_controls_; + + // TODO: we should only maintain one current control either from the well_state or from well_controls struct. + // Either one can be more favored depending on the final strategy for the initilzation of the well control + const int current = well_state.currentControls()[index_of_well_]; + const double target_rate = well_controls_iget_target(wc, current); + + std::vector g = {1,1,0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + const double* distr = well_controls_iget_distr(wc, current); + for (int p = 0; p < np; ++p) { + if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value + F[p] /= distr[p]; + } else { + F[p] = 0.; + } + } + } else { + for (int p = 0; p < np; ++p) { + F[p] /= g[p]; + } + } + + switch (well_controls_iget_type(wc, current)) { + case THP: // The BHP and THP both uses the total rate as first well variable. + case BHP: + { + primary_variables_[XvarWell] = xvar_well_old[XvarWell] - dwells[0][XvarWell]; + + switch (well_type_) { + case INJECTOR: + for (int p = 0; p < np; ++p) { + const double comp_frac = comp_frac_[p]; + well_state.wellRates()[index_of_well_ * np + p] = comp_frac * primary_variables_[XvarWell]; + } + break; + case PRODUCER: + for (int p = 0; p < np; ++p) { + well_state.wellRates()[index_of_well_ * np + p] = primary_variables_[XvarWell] * F[p]; + } + break; + } + + if (well_controls_iget_type(wc, current) == THP) { + + // Calculate bhp from thp control and well rates + std::vector rates(3, 0.0); // the vfp related only supports three phases for the moment + + const Opm::PhaseUsage& pu = phaseUsage(); + if (active()[ Water ]) { + rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ]; + } + if (active()[ Oil ]) { + rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; + } + if (active()[ Gas ]) { + rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ]; + } + + well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current); + } + } + break; + case SURFACE_RATE: // Both rate controls use bhp as first well variable + case RESERVOIR_RATE: + { + const int sign1 = dwells[0][XvarWell] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[0][XvarWell]),std::abs(xvar_well_old[XvarWell])*dBHPLimit); + primary_variables_[XvarWell] = std::max(xvar_well_old[XvarWell] - dx1_limited,1e5); + well_state.bhp()[index_of_well_] = primary_variables_[XvarWell]; + + if (well_controls_iget_type(wc, current) == SURFACE_RATE) { + if (well_type_ == PRODUCER) { + + const double* distr = well_controls_iget_distr(wc, current); + + double F_target = 0.0; + for (int p = 0; p < np; ++p) { + F_target += distr[p] * F[p]; + } + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate / F_target; + } + } else { + + for (int p = 0; p < np; ++p) { + well_state.wellRates()[index_of_well_ * np + p] = comp_frac_[p] * target_rate; + } + } + } else { // RESERVOIR_RATE + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate; + } + } + } + break; + } // end of switch (well_controls_iget_type(wc, current)) + + // for the wells having a THP constaint, we should update their thp value + // If it is under THP control, it will be set to be the target value. Otherwise, + // the thp value will be calculated based on the bhp value, assuming the bhp value is correctly calculated. + const int nwc = well_controls_get_num(wc); + // Looping over all controls until we find a THP constraint + int ctrl_index = 0; + for ( ; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(wc, ctrl_index) == THP) { + // the current control + const int current = well_state.currentControls()[index_of_well_]; + // If under THP control at the moment + if (current == ctrl_index) { + const double thp_target = well_controls_iget_target(wc, current); + well_state.thp()[index_of_well_] = thp_target; + } else { // otherwise we calculate the thp from the bhp value + const Opm::PhaseUsage& pu = phaseUsage(); + std::vector rates(3, 0.0); + + if (active()[ Water ]) { + rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ]; + } + if (active()[ Oil ]) { + rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ]; + } + if (active()[ Gas ]) { + rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ]; + } + + const double bhp = well_state.bhp()[index_of_well_]; + + well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp); + } + + // the THP control is found, we leave the loop now + break; + } + } // end of for loop for seaching THP constraints + + // no THP constraint found + if (ctrl_index == nwc) { // not finding a THP contstraints + well_state.thp()[index_of_well_] = 0.0; + } + } + + + + + + template + void + StandardWell:: + updateWellStateWithTarget(const int current, + WellState& xw) const + { + // number of phases + const int np = number_of_phases_; + const int well_index = index_of_well_; + const WellControls* wc = well_controls_; + // Updating well state and primary variables. + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[well_index] = target; + // TODO: similar to the way below to handle THP + // we should not something related to thp here when there is thp constraint + break; + + case THP: { + xw.thp()[well_index] = target; + + const Opm::PhaseUsage& pu = phaseUsage(); + std::vector rates(3, 0.0); + if (active()[ Water ]) { + rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; + } + if (active()[ Oil ]) { + rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; + } + if (active()[ Gas ]) { + rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ]; + } + + const int table_id = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + xw.bhp()[well_index] = calculateBhpFromThp(rates, current); + break; + } + + case RESERVOIR_RATE: // intentional fall-through + case SURFACE_RATE: + // checking the number of the phases under control + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + + assert(numPhasesWithTargetsUnderThisControl > 0); + + if (well_type_ == INJECTOR) { + // assign target value as initial guess for injectors + // only handles single phase control at the moment + assert(numPhasesWithTargetsUnderThisControl == 1); + + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.) { + xw.wellRates()[np*well_index + phase] = target / distr[phase]; + } else { + xw.wellRates()[np * well_index + phase] = 0.; + } + } + } else if (well_type_ == PRODUCER) { + // update the rates of phases under control based on the target, + // and also update rates of phases not under control to keep the rate ratio, + // assuming the mobility ratio does not change for the production wells + double original_rates_under_phase_control = 0.0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase]; + } + } + + if (original_rates_under_phase_control != 0.0 ) { + double scaling_factor = target / original_rates_under_phase_control; + + for (int phase = 0; phase < np; ++phase) { + xw.wellRates()[np * well_index + phase] *= scaling_factor; + } + } else { // scaling factor is not well defied when original_rates_under_phase_control is zero + // separating targets equally between phases under control + const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase]; + } else { + // this only happens for SURFACE_RATE control + xw.wellRates()[np * well_index + phase] = target_rate_divided; + } + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + break; + } // end of switch + + updatePrimaryVariables(xw); + } + + + + + + template + void + StandardWell:: + updateWellControl(WellState& xw, + wellhelpers::WellSwitchingLogger& logger) const + { + const int np = number_of_phases_; + const int nw = xw.bhp().size(); + const int w = index_of_well_; + + const int old_control_index = xw.currentControls()[w]; + + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + WellControls* wc = well_controls_; + + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + // the current control index + int current = xw.currentControls()[w]; + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, well_type_, wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + well_controls_set_current( wc, current); + } + + // the new well control indices after all the related updates, + const int updated_control_index = xw.currentControls()[w]; + + // checking whether control changed + if (updated_control_index != old_control_index) { + logger.wellSwitched(name(), + well_controls_iget_type(wc, old_control_index), + well_controls_iget_type(wc, updated_control_index)); + } + + if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) { + updateWellStateWithTarget(updated_control_index, xw); + } + } + + + + + + template + void + StandardWell:: + 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 + { + const int nperf = number_of_perforations_; + // TODO: can make this a member? + const int nw = xw.bhp().size(); + const int numComp = numComponents(); + const PhaseUsage& pu = *phase_usage_; + b_perf.resize(nperf*numComp); + surf_dens_perf.resize(nperf*numComp); + const int w = index_of_well_; + + //rs and rv are only used if both oil and gas is present + if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + } + + // Compute the average pressure in each well block + for (int perf = 0; perf < nperf; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + // TODO: this is another place to show why WellState need to be a vector of WellState. + // TODO: to check why should be perf - 1 + const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1]; + const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2; + const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); + + if (pu.phase_used[BlackoilPhases::Aqua]) { + 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 * numComp; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + 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) { + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + 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 double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + if (gasrate > 0) { + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + // Surface density. + for (int p = 0; p < pu.num_phases; ++p) { + surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + } + + // We use cell values for solvent injector + if (has_solvent) { + b_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventRefDensity(); + } + } + } + + + + + + template + void + StandardWell:: + computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf) + { + // Verify that we have consistent input. + const int np = number_of_phases_; + const int nperf = number_of_perforations_; + const int num_comp = numComponents(); + const PhaseUsage* phase_usage = phase_usage_; + + // 1. Compute the flow (in surface volume units for each + // component) exiting up the wellbore from each perforation, + // taking into account flow from lower in the well, and + // in/out-flow at each perforation. + std::vector q_out_perf(nperf*num_comp); + + // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore + // Iterate over well perforations from bottom to top. + for (int perf = nperf - 1; perf >= 0; --perf) { + for (int component = 0; component < num_comp; ++component) { + if (perf == nperf - 1) { + // This is the bottom perforation. No flow from below. + q_out_perf[perf*num_comp+ component] = 0.0; + } else { + // Set equal to flow from below. + q_out_perf[perf*num_comp + component] = q_out_perf[(perf+1)*num_comp + component]; + } + // Subtract outflow through perforation. + q_out_perf[perf*num_comp + component] -= perfComponentRates[perf*num_comp + component]; + } + } + + // 2. Compute the component mix at each perforation as the + // absolute values of the surface rates divided by their sum. + // Then compute volume ratios (formation factors) for each perforation. + // Finally compute densities for the segments associated with each perforation. + const int gaspos = phase_usage->phase_pos[BlackoilPhases::Vapour]; + const int oilpos = phase_usage->phase_pos[BlackoilPhases::Liquid]; + std::vector mix(num_comp,0.0); + std::vector x(num_comp); + std::vector surf_dens(num_comp); + std::vector dens(nperf); + + for (int perf = 0; perf < nperf; ++perf) { + // Find component mix. + const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf, + q_out_perf.begin() + num_comp*(perf+1), 0.0); + if (tot_surf_rate != 0.0) { + for (int component = 0; component < num_comp; ++component) { + mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); + } + } else { + // No flow => use well specified fractions for mix. + for (int phase = 0; phase < np; ++phase) { + mix[phase] = comp_frac_[phase]; + } + // intialize 0.0 for comIdx >= np; + } + // Compute volume ratio. + x = mix; + double rs = 0.0; + double rv = 0.0; + if (!rsmax_perf.empty() && mix[oilpos] > 0.0) { + rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); + } + if (!rvmax_perf.empty() && mix[gaspos] > 0.0) { + rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); + } + if (rs != 0.0) { + // Subtract gas in oil from gas mixture + x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv); + } + if (rv != 0.0) { + // Subtract oil in gas from oil mixture + x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);; + } + double volrat = 0.0; + for (int component = 0; component < num_comp; ++component) { + volrat += x[component] / b_perf[perf*num_comp+ component]; + } + for (int component = 0; component < num_comp; ++component) { + surf_dens[component] = surf_dens_perf[perf*num_comp+ component]; + } + + // Compute segment density. + perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; + } + } + + + + + template + void + StandardWell:: + computeConnectionPressureDelta() + { + // Algorithm: + + // We'll assume the perforations are given in order from top to + // bottom for each well. By top and bottom we do not necessarily + // mean in a geometric sense (depth), but in a topological sense: + // the 'top' perforation is nearest to the surface topologically. + // Our goal is to compute a pressure delta for each perforation. + + // 1. Compute pressure differences between perforations. + // dp_perf will contain the pressure difference between a + // perforation and the one above it, except for the first + // perforation for each well, for which it will be the + // difference to the reference (bhp) depth. + + const int nperf = number_of_perforations_; + perf_pressure_diffs_.resize(nperf, 0.0); + + for (int perf = 0; perf < nperf; ++perf) { + const double z_above = perf == 0 ? ref_depth_ : perf_depth_[perf - 1]; + const double dz = perf_depth_[perf] - z_above; + perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * gravity_; + } + + // 2. Compute pressure differences to the reference point (bhp) by + // accumulating the already computed adjacent pressure + // differences, storing the result in dp_perf. + // This accumulation must be done per well. + const auto beg = perf_pressure_diffs_.begin(); + const auto end = perf_pressure_diffs_.end(); + std::partial_sum(beg, end, beg); + } + + + + + + template + typename StandardWell::ConvergenceReport + StandardWell:: + getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg, + const ModelParameters& param) const + { + typedef double Scalar; + typedef std::vector< Scalar > Vector; + + const int np = number_of_phases_; + const int numComp = numComponents(); + + // the following implementation assume that the polymer is always after the w-o-g phases + // For the polymer case, there is one more mass balance equations of reservoir than wells + assert((int(B_avg.size()) == numComp) || has_polymer); + + const double tol_wells = param.tolerance_wells_; + const double maxResidualAllowed = param.max_residual_allowed_; + + // TODO: it should be the number of numWellEq + // using numComp here for flow_ebos running 2p case. + std::vector res(numComp); + for (int comp = 0; comp < numComp; ++comp) { + // magnitude of the residual matters + res[comp] = std::abs(resWell_[0][comp]); + } + + Vector well_flux_residual(numComp); + + // Finish computation + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + { + well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx]; + } + + ConvergenceReport report; + // checking if any NaN or too large residuals found + // TODO: not understand why phase here while component in other places. + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + + if (std::isnan(well_flux_residual[phaseIdx])) { + report.nan_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; + report.nan_residual_wells.push_back(problem_well); + } else { + if (well_flux_residual[phaseIdx] > maxResidualAllowed) { + report.too_large_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; + report.too_large_residual_wells.push_back(problem_well); + } + } + } + + if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found + // check convergence + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + { + report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells); + } + } else { // abnormal values found and no need to check the convergence + report.converged = false; + } + + return report; + } + + + + + + template + void + StandardWell:: + computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf) + { + // Compute densities + const int nperf = number_of_perforations_; + const int numComponent = numComponents(); + const int np = number_of_phases_; + std::vector perfRates(b_perf.size(),0.0); + + for (int perf = 0; perf < nperf; ++perf) { + for (int phase = 0; phase < np; ++phase) { + perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase]; + } + if(has_solvent) { + perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf]; + } + } + + computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + computeConnectionPressureDelta(); + + } + + + + + + template + void + StandardWell:: + computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& well_state) + { + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + computeWellConnectionDensitesPressures(well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + } + + + + + + template + void + StandardWell:: + solveEqAndUpdateWellState(const ModelParameters& param, + WellState& well_state) + { + // We assemble the well equations, then we check the convergence, + // which is why we do not put the assembleWellEq here. + BVectorWell dx_well(1); + invDuneD_.mv(resWell_, dx_well); + + updateWellState(dx_well, param, well_state); + } + + + + + + template + void + StandardWell:: + computeAccumWell() + { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value(); + } + } + + + + + + template + void + StandardWell:: + apply(const BVector& x, BVector& Ax) const + { + assert( Bx_.size() == duneB_.N() ); + assert( invDrw_.size() == invDuneD_.N() ); + + // Bx_ = duneB_ * x + duneB_.mv(x, Bx_); + // invDBx = invDuneD_ * Bx_ + // TODO: with this, we modified the content of the invDrw_. + // Is it necessary to do this to save some memory? + BVectorWell& invDBx = invDrw_; + invDuneD_.mv(Bx_, invDBx); + + // Ax = Ax - duneC_^T * invDBx + duneC_.mmtv(invDBx,Ax); + } + + + + + template + void + StandardWell:: + apply(BVector& r) const + { + assert( invDrw_.size() == invDuneD_.N() ); + + // invDrw_ = invDuneD_ * resWell_ + invDuneD_.mv(resWell_, invDrw_); + // r = r - duneC_^T * invDrw_ + duneC_.mmtv(invDrw_, r); + } + + + + + + template + void + StandardWell:: + recoverSolutionWell(const BVector& x, BVectorWell& xw) const + { + BVectorWell resWell = resWell_; + // resWell = resWell - B * x + duneB_.mmv(x, resWell); + // xw = D^-1 * resWell + invDuneD_.mv(resWell, xw); + } + + + + + + template + void + StandardWell:: + recoverWellSolutionAndUpdateWellState(const BVector& x, + const ModelParameters& param, + WellState& well_state) const + { + BVectorWell xw(1); + recoverSolutionWell(x, xw); + updateWellState(xw, param, well_state); + } + + + + + + template + void + StandardWell:: + computeWellRatesWithBhp(const Simulator& ebosSimulator, + const EvalWell& bhp, + std::vector& well_flux) const + { + const int np = number_of_phases_; + const int numComp = numComponents(); + well_flux.resize(np, 0.0); + + const bool allow_cf = crossFlowAllowed(ebosSimulator); + + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + // flux for each perforation + std::vector cq_s(numComp, 0.0); + std::vector mob(numComp, 0.0); + getMobility(ebosSimulator, perf, mob); + computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + + for(int p = 0; p < np; ++p) { + well_flux[p] += cq_s[p].value(); + } + } + } + + + + + + template + std::vector + StandardWell:: + computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const double initial_bhp, // bhp from BHP constraints + const std::vector& initial_potential) const + { + // TODO: pay attention to the situation that finally the potential is calculated based on the bhp control + // TODO: should we consider the bhp constraints during the iterative process? + const int np = number_of_phases_; + + assert( np == int(initial_potential.size()) ); + + std::vector potentials = initial_potential; + std::vector old_potentials = potentials; // keeping track of the old potentials + + double bhp = initial_bhp; + double old_bhp = bhp; + + bool converged = false; + const int max_iteration = 1000; + const double bhp_tolerance = 1000.; // 1000 pascal + + int iteration = 0; + + while ( !converged && iteration < max_iteration ) { + // for each iteration, we calculate the bhp based on the rates/potentials with thp constraints + // with considering the bhp value from the bhp limits. At the beginning of each iteration, + // we initialize the bhp to be the bhp value from the bhp limits. Then based on the bhp values calculated + // from the thp constraints, we decide the effective bhp value for well potential calculation. + bhp = initial_bhp; + + // The number of the well controls/constraints + const int nwc = well_controls_get_num(well_controls_); + + for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(well_controls_, ctrl_index) == THP) { + const Opm::PhaseUsage& pu = *phase_usage_; + + std::vector rates(3, 0.0); + if (active()[ Water ]) { + rates[ Water ] = potentials[pu.phase_pos[ Water ] ]; + } + if (active()[ Oil ]) { + rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ]; + } + if (active()[ Gas ]) { + rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ]; + } + + const double bhp_calculated = calculateBhpFromThp(rates, ctrl_index); + + if (well_type_ == INJECTOR && bhp_calculated < bhp ) { + bhp = bhp_calculated; + } + + if (well_type_ == PRODUCER && bhp_calculated > bhp) { + bhp = bhp_calculated; + } + } + } + + // there should be always some available bhp/thp constraints there + if (std::isinf(bhp) || std::isnan(bhp)) { + OPM_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << name()); + } + + converged = std::abs(old_bhp - bhp) < bhp_tolerance; + + computeWellRatesWithBhp(ebosSimulator, bhp, potentials); + + // checking whether the potentials have valid values + for (const double value : potentials) { + if (std::isinf(value) || std::isnan(value)) { + OPM_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << name()); + } + } + + if (!converged) { + old_bhp = bhp; + for (int p = 0; p < np; ++p) { + // TODO: improve the interpolation, will it always be valid with the way below? + // TODO: finding better paramters, better iteration strategy for better convergence rate. + const double potential_update_damping_factor = 0.001; + potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p]; + old_potentials[p] = potentials[p]; + } + } + + ++iteration; + } + + if (!converged) { + OPM_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << name()); + } + + return potentials; + } + + + + + + template + void + StandardWell:: + computeWellPotentials(const Simulator& ebosSimulator, + const WellState& well_state, + std::vector& well_potentials) const + { + const int np = number_of_phases_; + + well_potentials.resize(np, 0.0); + + // get the bhp value based on the bhp constraints + const double bhp = mostStrictBhpFromBhpLimits(); + + // does the well have a THP related constraint? + if ( !wellHasTHPConstraints() ) { + assert(std::abs(bhp) != std::numeric_limits::max()); + + computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials); + } else { + // the well has a THP related constraint + // checking whether a well is newly added, it only happens at the beginning of the report step + if ( !well_state.isNewWell(index_of_well_) ) { + for (int p = 0; p < np; ++p) { + // This is dangerous for new added well + // since we are not handling the initialization correctly for now + well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p]; + } + } else { + // We need to generate a reasonable rates to start the iteration process + computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials); + for (double& value : well_potentials) { + // make the value a little safer in case the BHP limits are default ones + // TODO: a better way should be a better rescaling based on the investigation of the VFP table. + const double rate_safety_scaling_factor = 0.00001; + value *= rate_safety_scaling_factor; + } + } + + well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials); + } + } + + + + + + template + void + StandardWell:: + updatePrimaryVariables(const WellState& well_state) const + { + const int np = number_of_phases_; + const int well_index = index_of_well_; + const WellControls* wc = well_controls_; + const double* distr = well_controls_get_current_distr(wc); + + std::vector g = {1.0, 1.0, 0.01}; + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + for (int phase = 0; phase < np; ++phase) { + g[phase] = distr[phase]; + } + } + + switch (well_controls_get_current_type(wc)) { + case THP: + case BHP: { + primary_variables_[XvarWell] = 0.0; + if (well_type_ == INJECTOR) { + for (int p = 0; p < np; ++p) { + primary_variables_[XvarWell] += well_state.wellRates()[np*well_index + p] * comp_frac_[p]; + } + } else { + for (int p = 0; p < np; ++p) { + primary_variables_[XvarWell] += g[p] * well_state.wellRates()[np*well_index + p]; + } + } + break; + } + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + primary_variables_[XvarWell] = well_state.bhp()[well_index]; + break; + } // end of switch + + double tot_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + tot_well_rate += g[p] * well_state.wellRates()[np*well_index + p]; + } + if(std::abs(tot_well_rate) > 0) { + if (active()[ Water ]) { + primary_variables_[WFrac] = g[Water] * well_state.wellRates()[np*well_index + Water] / tot_well_rate; + } + if (active()[ Gas ]) { + primary_variables_[GFrac] = g[Gas] * (well_state.wellRates()[np*well_index + Gas] - well_state.solventWellRate(well_index)) / tot_well_rate ; + } + if (has_solvent) { + primary_variables_[SFrac] = g[Gas] * well_state.solventWellRate(well_index) / tot_well_rate ; + } + } else { // tot_well_rate == 0 + if (well_type_ == INJECTOR) { + // only single phase injection handled + if (active()[Water]) { + if (distr[Water] > 0.0) { + primary_variables_[WFrac] = 1.0; + } else { + primary_variables_[WFrac] = 0.0; + } + } + + if (active()[Gas]) { + if (distr[Gas] > 0.0) { + primary_variables_[GFrac] = 1.0 - wsolvent(); + if (has_solvent) { + primary_variables_[SFrac] = wsolvent(); + } + } else { + primary_variables_[GFrac] = 0.0; + } + } + + // TODO: it is possible to leave injector as a oil well, + // when F_w and F_g both equals to zero, not sure under what kind of circumstance + // this will happen. + } else if (well_type_ == PRODUCER) { // producers + // TODO: the following are not addressed for the solvent case yet + if (active()[Water]) { + primary_variables_[WFrac] = 1.0 / np; + } + if (active()[Gas]) { + primary_variables_[GFrac] = 1.0 / np; + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + } + } + + + + + + template + template + ValueType + StandardWell:: + calculateBhpFromThp(const std::vector& rates, + const int control_index) const + { + // TODO: when well is under THP control, the BHP is dependent on the rates, + // the well rates is also dependent on the BHP, so it might need to do some iteration. + // However, when group control is involved, change of the rates might impacts other wells + // so iterations on a higher level will be required. Some investigation might be needed when + // we face problems under THP control. + + assert(int(rates.size()) == 3); // the vfp related only supports three phases now. + + const ValueType aqua = rates[Water]; + const ValueType liquid = rates[Oil]; + const ValueType vapour = rates[Gas]; + + const int vfp = well_controls_iget_vfp(well_controls_, control_index); + const double& thp = well_controls_iget_target(well_controls_, control_index); + const double& alq = well_controls_iget_alq(well_controls_, control_index); + + // pick the density in the top layer + const double rho = perf_densities_[0]; + + ValueType bhp = 0.; + if (well_type_ == INJECTOR) { + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(); + + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); + + bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type_ == PRODUCER) { + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(); + + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); + + bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + return bhp; + } + + + + + + template + double + StandardWell:: + calculateThpFromBhp(const std::vector& rates, + const int control_index, + const double bhp) const + { + assert(int(rates.size()) == 3); // the vfp related only supports three phases now. + + const double aqua = rates[Water]; + const double liquid = rates[Oil]; + const double vapour = rates[Gas]; + + const int vfp = well_controls_iget_vfp(well_controls_, control_index); + const double& alq = well_controls_iget_alq(well_controls_, control_index); + + // pick the density in the top layer + const double rho = perf_densities_[0]; + + double thp = 0.0; + if (well_type_ == INJECTOR) { + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(); + + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); + + thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp); + } + else if (well_type_ == PRODUCER) { + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(); + + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); + + thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + return thp; + } + +} diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 58674c1db..48114afa5 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -39,8 +39,6 @@ #include #include #include -#include -#include #include #include #include @@ -49,27 +47,19 @@ #include #include #include +#include +#include #include #include #include #include -#include #include -#include namespace Opm { -enum WellVariablePositions { - XvarWell = 0, - WFrac = 1, - GFrac = 2, - SFrac = 3 -}; - - /// Class for handling the standard well model. template class StandardWellsDense { @@ -82,26 +72,17 @@ enum WellVariablePositions { 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 typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; static const int numEq = BlackoilIndices::numEq; - static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // //numEq; //number of wellEq is only numEq for polymer - static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; + // TODO: where we should put these types, WellInterface or Well Model? + // or there is some other strategy, like TypeTag typedef Dune::FieldVector VectorBlockType; - typedef Dune::FieldMatrix MatrixBlockType; - typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; - typedef DenseAd::Evaluation EvalWell; - typedef DenseAd::Evaluation Eval; typedef Ewoms::BlackOilPolymerModule PolymerModule; // For the conversion between the surface volume rate and resrevoir voidage rate @@ -121,13 +102,143 @@ enum WellVariablePositions { const std::vector& active_arg, const double gravity_arg, const std::vector& depth_arg, - const std::vector& pv_arg, long int global_nc, const Grid& grid); void setVFPProperties(const VFPProperties* vfp_properties_arg); - /// The number of components in the model. + + SimulatorReport assemble(Simulator& ebosSimulator, + const int iterationIdx, + const double dt, + WellState& well_state); + + // substract Binv(D)rw from r; + void apply( BVector& r) const; + + // subtract B*inv(D)*C * x from A*x + void apply(const BVector& x, BVector& Ax) const; + + // apply well model with scaling of alpha + void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const; + + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const; + + int numWells() const; + + /// return true if wells are available in the reservoir + bool wellsActive() const; + + void setWellsActive(const bool wells_active); + + /// return true if wells are available on this process + bool localWellsActive() const; + + bool getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg) 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; + + WellCollection* wellCollection() const; + + + protected: + bool wells_active_; + const Wells* wells_; + const std::vector< const Well* > wells_ecl_; + + // the number of wells in this process + // trying not to use things from Wells struct + // TODO: maybe a better name to emphasize it is local? + const int number_of_wells_; + + const int number_of_phases_; + + using WellInterfacePtr = std::unique_ptr >; + // a vector of all the wells. + // eventually, the wells_ above should be gone. + // the name is just temporary + // later, might make share_ptr const later. + std::vector well_container_; + + using ConvergenceReport = typename WellInterface::ConvergenceReport; + + // create the well container + static std::vector createWellContainer(const Wells* wells, + const std::vector& wells_ecl, + const int time_step); + + // Well collection is used to enforce the group control + WellCollection* well_collection_; + + ModelParameters param_; + bool terminal_output_; + bool has_solvent_; + bool has_polymer_; + int current_timeIdx_; + + PhaseUsage phase_usage_; + std::vector active_; + const RateConverterType& rate_converter_; + + // the number of the cells in the local grid + int number_of_cells_; + + long int global_nc_; + + // used to better efficiency of calcuation + mutable BVector scaleAddRes_; + + void updateWellControls(WellState& xw) const; + + void updateGroupControls(WellState& well_state) const; + + // setting the well_solutions_ based on well_state. + void updatePrimaryVariables(const WellState& well_state) const; + + void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const; + + void computeRepRadiusPerfLength(const Grid& grid); + + + void computeAverageFormationFactor(Simulator& ebosSimulator, + std::vector& B_avg) const; + + void applyVREPGroupControl(WellState& well_state) const; + + void computeWellVoidageRates(const WellState& well_state, + std::vector& well_voidage_rates, + std::vector& voidage_conversion_coeffs) const; + + // Calculating well potentials for each well + // TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP. + void computeWellPotentials(const Simulator& ebosSimulator, + const WellState& well_state, + std::vector& well_potentials) const; + + const std::vector& wellPerfEfficiencyFactors() const; + + void calculateEfficiencyFactors(); + + void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw) const; + + SimulatorReport solveWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state) const; + + void computeAccumWells() const; + + void initPrimaryVariablesEvaluation() const; + + // The number of components in the model. int numComponents() const { if (numPhases() == 2) { @@ -141,286 +252,21 @@ enum WellVariablePositions { return numComp; } + int numPhases() const; - SimulatorReport assemble(Simulator& ebosSimulator, - const int iterationIdx, - const double dt, - WellState& well_state); + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + + void resetWellControlFromState(const WellState& xw) const; void assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, - bool only_wells); + bool only_wells) const; - void - getMobility(const Simulator& ebosSimulator, - const int w, - const int perf, - const int cell_idx, - std::vector& mob) const; - - bool allow_cross_flow(const int w, const Simulator& ebosSimulator) const; - - void localInvert(Mat& istlA) const; - - void print(Mat& istlA) const; - - // substract Binv(D)rw from r; - void apply( BVector& r) const; - - // subtract B*inv(D)*C * x from A*x - void apply(const BVector& x, BVector& Ax) const; - - // apply well model with scaling of alpha - void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const; - - // xw = inv(D)*(rw - C*x) - void recoverVariable(const BVector& x, BVector& xw) const; - - int flowPhaseToEbosCompIdx( const int phaseIdx ) const; - - int flowToEbosPvIdx( const int flowPv ) const; - - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - - std::vector - extractPerfData(const std::vector& in) const; - - int numPhases() const; - - int numCells() const; - - void resetWellControlFromState(const WellState& xw) const; - - const Wells& wells() const; - - const Wells* wellsPointer() const; - - /// return true if wells are available in the reservoir - bool wellsActive() const; - - void setWellsActive(const bool wells_active); - - /// return true if wells are available on this process - bool localWellsActive() const; - - int numWellVars() const; - - /// Density of each well perforation - const std::vector& wellPerforationDensities() const; - - /// Diff to bhp for each well perforation. - const std::vector& wellPerforationPressureDiffs() const; - - EvalWell extendEval(const Eval& in) const; - - void setWellVariables(const WellState& xw); - - void print(const EvalWell& in) const; - - void computeAccumWells(); - - 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; - - SimulatorReport solveWellEq(Simulator& ebosSimulator, - const double dt, - WellState& well_state); - - void printIf(const int c, const double x, const double y, const double eps, const std::string type) const; - - std::vector residual() const; - - bool getWellConvergence(Simulator& ebosSimulator, - const int iteration) const; - - void computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw); - - 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; - - - - 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 computeWellConnectionDensitesPressures(const WellState& xw, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf, - const std::vector& depth_perf, - const double grav); - - - // Calculating well potentials for each well - // TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP. - void computeWellPotentials(const Simulator& ebosSimulator, - const WellState& well_state, - std::vector& well_potentials) const; - - // TODO: some preparation work, mostly related to group control and RESV, + // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) void prepareTimeStep(const Simulator& ebos_simulator, WellState& well_state); - - WellCollection* wellCollection() const; - - const std::vector& - wellPerfEfficiencyFactors() const; - - void calculateEfficiencyFactors(); - - void computeWellVoidageRates(const WellState& well_state, - std::vector& well_voidage_rates, - std::vector& voidage_conversion_coeffs) const; - - void applyVREPGroupControl(WellState& well_state) const; - - protected: - bool wells_active_; - const Wells* wells_; - const std::vector< const Well* > wells_ecl_; - - // Well collection is used to enforce the group control - WellCollection* well_collection_; - - ModelParameters param_; - bool terminal_output_; - bool has_solvent_; - bool has_polymer_; - int current_timeIdx_; - - PhaseUsage phase_usage_; - std::vector active_; - const VFPProperties* vfp_properties_; - double gravity_; - const RateConverterType& rate_converter_; - - // The efficiency factor for each connection. It is specified based on wells and groups, - // We calculate the factor for each connection for the computation of contributions to the mass balance equations. - // By default, they should all be one. - std::vector well_perforation_efficiency_factors_; - // the depth of the all the cell centers - // for standard Wells, it the same with the perforation depth - std::vector cell_depths_; - std::vector pv_; - - std::vector well_perforation_densities_; - std::vector well_perforation_pressure_diffs_; - - std::vector wpolymer_; - std::vector wsolvent_; - - std::vector wells_rep_radius_; - std::vector wells_perf_length_; - std::vector wells_bore_diameter_; - - std::vector wellVariables_; - std::vector F0_; - - Mat duneB_; - Mat duneC_; - Mat invDuneD_; - - BVector resWell_; - - long int global_nc_; - - mutable BVector Cx_; - mutable BVector invDrw_; - mutable BVector scaleAddRes_; - - double dbhpMaxRel() const {return param_.dbhp_max_rel_; } - double dWellFractionMax() const {return param_.dwell_fraction_max_; } - - // protected methods - EvalWell getBhp(const int wellIdx) const; - - EvalWell getQs(const int wellIdx, const int compIdx) const; - - EvalWell wellVolumeFraction(const int wellIdx, const int compIdx) 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 compIdx) const; - - bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const int well_number) const; - - using WellMapType = typename WellState::WellMapType; - using WellMapEntryType = typename WellState::mapentry_t; - - // a tuple type for ratio limit check. - // first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three - // values should not be used. - // second value indicates whehter there is only one connection left. - // third value indicates the indx of the worst-offending connection. - // the last value indicates the extent of the violation for the worst-offending connection, which is defined by - // the ratio of the actual value to the value of the violated limit. - using RatioCheckTuple = std::tuple; - - enum ConnectionIndex { - INVALIDCONNECTION = -10000 - }; - - - RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const; - - RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const; - - void updateWellStateWithTarget(const WellControls* wc, - const int current, - const int well_index, - WellState& xw) const; - - bool wellHasTHPConstraints(const int well_index) const; - - // TODO: maybe we should provide a light version of computeWellFlux, which does not include the - // calculation of the derivatives - 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 - std::vector - computeWellPotentialWithTHP(const Simulator& ebosSimulator, - const int well_index, - const double initial_bhp, // bhp from BHP constraints - const std::vector& initial_potential) const; - - double wsolvent(const int well_index) const; - - double wpolymer(const int well_index) const; - - void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const; - - void computeRepRadiusPerfLength(const Grid& grid); - - }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 11df056c3..ba841d2a5 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -15,6 +15,9 @@ namespace Opm { : wells_active_(wells_arg!=nullptr) , wells_(wells_arg) , wells_ecl_(wells_ecl) + , number_of_wells_(wells_arg ? (wells_arg->number_of_wells) : 0) + , number_of_phases_(wells_arg ? (wells_arg->number_of_phases) : 0) // TODO: not sure if it is proper for this way + , well_container_(createWellContainer(wells_arg, wells_ecl, current_timeIdx) ) , well_collection_(well_collection) , param_(param) , terminal_output_(terminal_output) @@ -22,18 +25,7 @@ namespace Opm { , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) , current_timeIdx_(current_timeIdx) , rate_converter_(rate_converter) - , 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 * numWellEq) : 0) - , F0_(wells_ ? (wells_arg->number_of_wells * numWellEq) : 0 ) { - if( wells_ ) - { - invDuneD_.setBuildMode( Mat::row_wise ); - duneC_.setBuildMode( Mat::row_wise ); - duneB_.setBuildMode( Mat::row_wise ); - } } @@ -47,7 +39,6 @@ namespace Opm { const std::vector& active_arg, const double gravity_arg, const std::vector& depth_arg, - const std::vector& pv_arg, long int global_nc, const Grid& grid) { @@ -60,22 +51,11 @@ namespace Opm { phase_usage_ = phase_usage_arg; active_ = active_arg; - gravity_ = gravity_arg; - cell_depths_ = extractPerfData(depth_arg); - pv_ = pv_arg; calculateEfficiencyFactors(); - // setup sparsity pattern for the matrices - //[A B^T [x = [ res - // C D] x_well] res_well] - - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const int nc = numCells(); - #ifndef NDEBUG - const auto pu = phase_usage_; + const auto& pu = phase_usage_; const int np = pu.num_phases; // assumes the gas fractions are stored after water fractions @@ -83,42 +63,6 @@ namespace Opm { assert (np == 3 || (np == 2 && !pu.phase_used[Gas]) ); #endif - // set invDuneD - invDuneD_.setSize( nw, nw, nw ); - - // set duneC - duneC_.setSize( nw, nc, nperf ); - - // set duneB - duneB_.setSize( nw, nc, nperf ); - - for (auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) { - // Add nonzeros for diagonal - row.insert(row.index()); - } - - for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) { - // Add nonzeros for diagonal - for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - row.insert(cell_idx); - } - } - - // make the B^T matrix - for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) { - for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - row.insert(cell_idx); - } - } - - resWell_.resize( nw ); - - // resize temporary class variables - Cx_.resize( duneC_.N() ); - invDrw_.resize( invDuneD_.N() ); - if (has_polymer_) { if (PolymerModule::hasPlyshlog()) { @@ -126,6 +70,13 @@ namespace Opm { } } + number_of_cells_ = Opm::UgGridHelpers::numCells(grid); + // do the initialization for all the wells + // TODO: to see whether we can postpone of the intialization of the well containers to + // optimize the usage of the following several member variables + for (auto& well : well_container_) { + well->init(&phase_usage_, &active_, depth_arg, gravity_arg, number_of_cells_); + } } @@ -137,7 +88,71 @@ namespace Opm { StandardWellsDense:: setVFPProperties(const VFPProperties* vfp_properties_arg) { - vfp_properties_ = vfp_properties_arg; + for (auto& well : well_container_) { + well->setVFPProperties(vfp_properties_arg); + } + } + + + + + + + template + int + StandardWellsDense:: + numWells() const + { + return number_of_wells_; + } + + + + + + template + std::vector::WellInterfacePtr > + StandardWellsDense:: + createWellContainer(const Wells* wells, + const std::vector< const Well* >& wells_ecl, + const int time_step) + { + std::vector well_container; + + const int nw = wells ? (wells->number_of_wells) : 0; + + if (nw > 0) { + well_container.reserve(nw); + + // With the following way, it will have the same order with wells struct + // Hopefully, it can generate the same residual history with master branch + for (int w = 0; w < nw; ++w) { + const std::string well_name = std::string(wells->name[w]); + + // finding the location of the well in wells_ecl + const int nw_wells_ecl = wells_ecl.size(); + int index_well = 0; + for (; index_well < nw_wells_ecl; ++index_well) { + if (well_name == wells_ecl[index_well]->name()) { + break; + } + } + + // It should be able to find in wells_ecl. + if (index_well == nw_wells_ecl) { + OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); + } + + const Well* well_ecl = wells_ecl[index_well]; + if (well_ecl->isMultiSegment(time_step)) { + OPM_THROW(Opm::NumericalProblem, "Not handling Multisegment Wells for now"); + } + + // Basically, we are handling all the wells as StandardWell for the moment + well_container.emplace_back(new StandardWell(well_ecl, time_step, wells) ); + } + } + return well_container; } @@ -163,8 +178,8 @@ namespace Opm { } updateWellControls(well_state); - // Set the primary variables for the wells - setWellVariables(well_state); + // Set the well primary variables based on the value of well solutions + initPrimaryVariablesEvaluation(); if (iterationIdx == 0) { computeWellConnectionPressures(ebosSimulator, well_state); @@ -191,278 +206,10 @@ namespace Opm { assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, - bool only_wells) + bool only_wells) const { - const int nw = wells().number_of_wells; - const int numComp = numComponents(); - const int np = numPhases(); - - // clear all entries - duneB_ = 0.0; - duneC_ = 0.0; - invDuneD_ = 0.0; - resWell_ = 0.0; - - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); - - const double volume = 0.002831684659200; // 0.1 cu ft; - for (int w = 0; w < nw; ++w) { - bool allow_cf = allow_cross_flow(w, ebosSimulator); - const EvalWell& bhp = getBhp(w); - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - - const int cell_idx = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - std::vector cq_s(numComp,0.0); - std::vector mob(numComp, 0.0); - getMobility(ebosSimulator, w, perf, cell_idx, mob); - computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - - 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[componentIdx] * well_perforation_efficiency_factors_[perf]; - - if (!only_wells) { - // subtract sum of component fluxes in the reservoir equation. - // need to consider the efficiency factor - ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value(); - } - - // subtract sum of phase fluxes in the well equations. - resWell_[w][componentIdx] -= cq_s[componentIdx].value(); - - // assemble the jacobians - for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix - } - invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); - } - - for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); - duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); - } - } - - // add a trivial equation for the dummy phase for 2p cases (Only support water + oil) - if ( numComp < numWellEq ) { - assert(!active_[ Gas ]); - invDuneD_[w][w][Gas][Gas] = 1.0; - } - - // Store the perforation phase flux for later usage. - if (has_solvent_ && componentIdx == solventSaturationIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) - well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value(); - } else { - well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value(); - } - } - - if (has_polymer_) { - EvalWell cq_s_poly = cq_s[Water]; - if (wells().type[w] == INJECTOR) { - cq_s_poly *= wpolymer(w); - } else { - cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); - } - if (!only_wells) { - for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx); - } - ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); - } - } - - // Store the perforation pressure for later usage. - well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; - } - - // add vol * dF/dt + Q to the well equations; - 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][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq); - } - resWell_[w][componentIdx] += resWell_loc.value(); - } - - // add trivial equation for polymer - if (has_polymer_) { - invDuneD_[w][w][contiPolymerEqIdx][polymerConcentrationIdx] = 1.0; // - } - } - - - - // do the local inversion of D. - localInvert( invDuneD_ ); - } - - - template - void - StandardWellsDense:: - getMobility(const Simulator& ebosSimulator, const int w, const int perf, const int cell_idx, std::vector& mob) const - { - - const int np = wells().number_of_phases; - assert (int(mob.size()) == numComponents()); - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - - // either use mobility of the perforation cell or calcualte its own - // based on passing the saturation table index - const int satid = wells().sat_table_id[perf] - 1; - const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell - - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); - } - if (has_solvent_) { - mob[solventSaturationIdx] = extendEval(intQuants.solventMobility()); - } - } else { - - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - Eval relativePerms[3] = { 0.0, 0.0, 0.0 }; - MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); - - // reset the satnumvalue back to original - materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); - - // compute the mobility - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); - } - - // this may not work if viscosity and relperms has been modified? - if (has_solvent_) { - OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent"); - } - } - - // modify the water mobility if polymer is present - if (has_polymer_) { - // assume fully mixture for wells. - EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); - - if (wells().type[w] == INJECTOR) { - const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex()); - mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) ); - } - - if (PolymerModule::hasPlyshlog()) { - // compute the well water velocity with out shear effects. - const int numComp = numComponents(); - bool allow_cf = allow_cross_flow(w, ebosSimulator); - const EvalWell& bhp = getBhp(w); - std::vector cq_s(numComp,0.0); - computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf]; - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - const auto& scaledDrainageInfo = - materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx); - const Scalar& Swcr = scaledDrainageInfo.Swcr; - const EvalWell poro = extendEval(intQuants.porosity()); - const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); - // guard against zero porosity and no water - const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); - EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); - - if (PolymerModule::hasShrate()) { - // TODO Use the same conversion as for the reservoar equations. - // Need the "permeability" of the well? - // For now use the same formula as in legacy. - waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf]; - } - EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); - EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, - intQuants.pvtRegionIndex(), - waterVelocity); - - // modify the mobility with the shear factor and recompute the well fluxes. - mob[ Water ] /= shearFactor; - } - } - - } - - - - - template - bool - StandardWellsDense:: - allow_cross_flow(const int w, const Simulator& ebosSimulator) const - { - if (wells().allow_cf[w]) { - return true; - } - - // check for special case where all perforations have cross flow - // then the wells must allow for cross flow - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - const int cell_idx = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); - EvalWell bhp = getBhp(w); - - // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + wellPerforationPressureDiffs()[perf]; - EvalWell drawdown = pressure - well_pressure; - - if (drawdown.value() < 0 && wells().type[w] == INJECTOR) { - return false; - } - - if (drawdown.value() > 0 && wells().type[w] == PRODUCER) { - return false; - } - } - return true; - } - - - - - - template - void - StandardWellsDense:: - localInvert(Mat& istlA) const - { - for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { - for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { - //std::cout << (*col) << std::endl; - (*col).invert(); - } - } - } - - - - - - template - void - StandardWellsDense:: - print(Mat& istlA) const - { - for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { - for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { - std::cout << row.index() << " " << col.index() << "/n \n"<<(*col) << std::endl; - } + for (int w = 0; w < number_of_wells_; ++w) { + well_container_[w]->assembleWellEq(ebosSimulator, dt, well_state, only_wells); } } @@ -470,6 +217,8 @@ namespace Opm { + // applying the well residual to reservoir residuals + // r = r - duneC_^T * invDuneD_ * resWell_ template void StandardWellsDense:: @@ -479,39 +228,36 @@ namespace Opm { return; } - assert( invDrw_.size() == invDuneD_.N() ); - - invDuneD_.mv(resWell_,invDrw_); - duneB_.mmtv(invDrw_, r); + for (auto& well : well_container_) { + well->apply(r); + } } + // Ax = A x - C D^-1 B x template void StandardWellsDense:: apply(const BVector& x, BVector& Ax) const { + // TODO: do we still need localWellsActive()? if ( ! localWellsActive() ) { return; } - assert( Cx_.size() == duneC_.N() ); - - BVector& invDCx = invDrw_; - assert( invDCx.size() == invDuneD_.N()); - - duneC_.mv(x, Cx_); - invDuneD_.mv(Cx_, invDCx); - duneB_.mmtv(invDCx,Ax); + for (auto& well : well_container_) { + well->apply(x, Ax); + } } + // Ax = Ax - alpha * C D^-1 B x template void StandardWellsDense:: @@ -526,7 +272,9 @@ namespace Opm { } scaleAddRes_ = 0.0; + // scaleAddRes_ = - C D^-1 B x apply( x, scaleAddRes_ ); + // Ax = Ax + alpha * scaleAddRes_ Ax.axpy( alpha, scaleAddRes_ ); } @@ -537,45 +285,11 @@ namespace Opm { template void StandardWellsDense:: - recoverVariable(const BVector& x, BVector& xw) const + recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const { - if ( ! localWellsActive() ) { - return; + for (auto& well : well_container_) { + well->recoverWellSolutionAndUpdateWellState(x, param_, well_state); } - BVector resWell = resWell_; - duneC_.mmv(x, resWell); - invDuneD_.mv(resWell, xw); - } - - - - - template - int - StandardWellsDense:: - flowToEbosPvIdx( const int flowPv ) const - { - const int flowToEbos[ 3 ] = { - BlackoilIndices::pressureSwitchIdx, - BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx - }; - - if (flowPv > 2 ) - return flowPv; - - return flowToEbos[ flowPv ]; - } - - template - int - StandardWellsDense:: - flowPhaseToEbosCompIdx( const int phaseIdx ) const - { - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; - if (phaseIdx > 2 ) - return phaseIdx; - return phaseToComp[ phaseIdx ]; } @@ -592,23 +306,6 @@ namespace Opm { return flowToEbos[ phaseIdx ]; } - template - std::vector - StandardWellsDense:: - extractPerfData(const std::vector& in) const - { - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - std::vector out(nperf); - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - const int well_idx = wells().well_cells[perf]; - out[perf] = in[well_idx]; - } - } - return out; - } - @@ -618,19 +315,7 @@ namespace Opm { StandardWellsDense:: numPhases() const { - return wells().number_of_phases; - } - - - - - - template - int - StandardWellsDense:: - numCells() const - { - return pv_.size(); + return number_of_phases_; } @@ -653,31 +338,6 @@ namespace Opm { - template - const Wells& - StandardWellsDense:: - wells() const - { - assert(wells_ != 0); - return *(wells_); - } - - - - - - template - const Wells* - StandardWellsDense:: - wellsPointer() const - { - return wells_; - } - - - - - template bool StandardWellsDense:: @@ -707,64 +367,7 @@ namespace Opm { StandardWellsDense:: localWellsActive() const { - return wells_ ? (wells_->number_of_wells > 0 ) : false; - } - - - - - - template - int - StandardWellsDense:: - numWellVars() const - { - if ( !localWellsActive() ) { - return 0; - } - - const int nw = wells().number_of_wells; - return numWellEq * nw; - } - - - - - - template - const std::vector& - StandardWellsDense:: - wellPerforationDensities() const - { - return well_perforation_densities_; - } - - - - - - template - const std::vector& - StandardWellsDense:: - wellPerforationPressureDiffs() const - { - return well_perforation_pressure_diffs_; - } - - - - - - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - extendEval(const Eval& in) const { - EvalWell out = 0.0; - out.setValue(in.value()); - for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { - out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx))); - } - return out; + return number_of_wells_ > 0; } @@ -774,22 +377,10 @@ namespace Opm { template void StandardWellsDense:: - setWellVariables(const WellState& xw) + initPrimaryVariablesEvaluation() const { - const int nw = wells().number_of_wells; - // for two-phase numComp < numWellEq - const int numComp = numComponents(); - for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) { - for (int w = 0; w < nw; ++w) { - const unsigned int idx = nw * eqIdx + w; - assert( idx < wellVariables_.size() ); - assert( idx < xw.wellSolutions().size() ); - EvalWell& eval = wellVariables_[ idx ]; - - eval = 0.0; - eval.setValue( xw.wellSolutions()[ idx ] ); - eval.setDerivative(numEq + eqIdx, 1.0); - } + for (auto& well : well_container_) { + well->initPrimaryVariablesEvaluation(); } } @@ -800,154 +391,10 @@ namespace Opm { template void StandardWellsDense:: - print(const EvalWell& in) const + computeAccumWells() const { - std::cout << in.value() << std::endl; - for (int i = 0; i < in.size; ++i) { - std::cout << in.derivative(i) << std::endl; - } - } - - - - - - template - void - StandardWellsDense:: - computeAccumWells() - { - const int nw = wells().number_of_wells; - for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { - for (int w = 0; w < nw; ++w) { - F0_[w + nw * eqIdx] = wellSurfaceVolumeFraction(w, eqIdx).value(); - } - } - } - - - - - - template - void - StandardWellsDense:: - 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 - { - const Opm::PhaseUsage& pu = phase_usage_; - const int np = wells().number_of_phases; - 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(numComp, 0.0); - for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); - b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); - } - if (has_solvent_) { - b_perfcells_dense[solventSaturationIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); - } - - // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + cdp; - EvalWell drawdown = pressure - well_pressure; - - // producing perforations - if ( drawdown.value() > 0 ) { - //Do nothing if crossflow is not allowed - if (!allow_cf && wells().type[w] == INJECTOR) { - return; - } - - // 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_sOil = cq_s[oilpos]; - const EvalWell cq_sGas = cq_s[gaspos]; - cq_s[gaspos] += rs * cq_sOil; - cq_s[oilpos] += rv * cq_sGas; - } - - } else { - //Do nothing if crossflow is not allowed - if (!allow_cf && wells().type[w] == PRODUCER) { - return; - } - - // Using total mobilities - EvalWell total_mob_dense = mob_perfcells_dense[0]; - for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) { - total_mob_dense += mob_perfcells_dense[componentIdx]; - } - - // injection perforations total volume rates - const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); - - // compute volume ratio between connection at standard conditions - EvalWell volumeRatio = 0.0; - if (active_[Water]) { - const int watpos = pu.phase_pos[Water]; - volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; - } - - if (has_solvent_) { - volumeRatio += cmix_s[solventSaturationIdx] / b_perfcells_dense[solventSaturationIdx]; - } - - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - - // Incorporate RS/RV factors if both oil and gas active - const EvalWell d = 1.0 - rv * rs; - - if (d.value() == 0.0) { - OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << wells().name[w] << " during flux calcuation" - << " with rs " << rs << " and rv " << rv); - } - - const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d; - //std::cout << "tmp_oil " <computeAccumWell(); } } @@ -960,16 +407,21 @@ namespace Opm { StandardWellsDense:: solveWellEq(Simulator& ebosSimulator, const double dt, - WellState& well_state) + WellState& well_state) const { - const int nw = wells().number_of_wells; + const int nw = number_of_wells_; WellState well_state0 = well_state; + const int numComp = numComponents(); + std::vector< Scalar > B_avg( numComp, Scalar() ); + computeAverageFormationFactor(ebosSimulator, B_avg); + int it = 0; bool converged; do { assembleWellEq(ebosSimulator, dt, well_state, true); - converged = getWellConvergence(ebosSimulator, it); + + converged = getWellConvergence(ebosSimulator, B_avg); // checking whether the group targets are converged if (wellCollection()->groupControlActive()) { @@ -983,10 +435,9 @@ namespace Opm { ++it; if( localWellsActive() ) { - BVector dx_well (nw); - invDuneD_.mv(resWell_, dx_well); - - updateWellState(dx_well, well_state); + for (auto& well : well_container_) { + well->solveEqAndUpdateWellState(param_, well_state); + } } // updateWellControls uses communication // Therefore the following is executed if there @@ -994,15 +445,20 @@ namespace Opm { if( wellsActive() ) { updateWellControls(well_state); - setWellVariables(well_state); + initPrimaryVariablesEvaluation(); } } while (it < 15); - if (!converged) { + if (converged) { + OpmLog::debug("Well equation solution gets converged with " + std::to_string(it) + " iterations"); + } else { + OpmLog::debug("Well equation solution failed in getting converged with " + std::to_string(it) + " iterations"); + well_state = well_state0; + updatePrimaryVariables(well_state); // also recover the old well controls for (int w = 0; w < nw; ++w) { - WellControls* wc = wells().ctrls[w]; + WellControls* wc = well_container_[w]->wellControls(); well_controls_set_current(wc, well_state.currentControls()[w]); } } @@ -1017,159 +473,59 @@ namespace Opm { - template - void - 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) { - std::cout << type << " " << c << ": "< - std::vector - StandardWellsDense:: - residual() const - { - if( ! wellsActive() ) - { - return std::vector(); - } - - const int nw = wells().number_of_wells; - const int numComp = numComponents(); - std::vector res(numEq*nw, 0.0); - for( int compIdx = 0; compIdx < numComp; ++compIdx) { - - for (int wellIdx = 0; wellIdx < nw; ++wellIdx) { - int idx = wellIdx + nw*compIdx; - res[idx] = resWell_[ wellIdx ][ compIdx ]; - } - } - return res; - } - - - - - template bool StandardWellsDense:: getWellConvergence(Simulator& ebosSimulator, - const int iteration) const + const std::vector& B_avg) const { - typedef double Scalar; - typedef std::vector< Scalar > Vector; + ConvergenceReport report; - 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( numComp, Scalar() ); - std::vector< Scalar > maxNormWell(numComp, Scalar() ); - - auto& grid = ebosSimulator.gridManager().grid(); - const auto& gridView = grid.leafGridView(); - ElementContext elemCtx(ebosSimulator); - const auto& elemEndIt = gridView.template end(); - - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; ++elemIt) - { - elemCtx.updatePrimaryStencil(*elemIt); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - - for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) - { - auto& B = B_avg[ phaseIdx ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); - - B += 1 / fs.invB(ebosPhaseIdx).value(); - } - if (has_solvent_) { - auto& B = B_avg[ solventSaturationIdx ]; - B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); - } + for (const auto& well : well_container_) { + report += well->getWellConvergence(ebosSimulator, B_avg, param_); } - // compute global average - grid.comm().sum(B_avg.data(), B_avg.size()); - for(auto& bval: B_avg) + // checking NaN residuals { - bval/=global_nc_; - } + bool nan_residual_found = report.nan_residual_found; + const auto& grid = ebosSimulator.gridManager().grid(); + int value = nan_residual_found ? 1 : 0; - auto res = residual(); - const int nw = res.size() / numComp; + nan_residual_found = grid.comm().max(value); - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) - { - for ( int w = 0; w < nw; ++w ) { - maxNormWell[compIdx] = std::max(maxNormWell[compIdx], std::abs(res[nw*compIdx + w])); - } - } - - grid.comm().max(maxNormWell.data(), maxNormWell.size()); - - Vector well_flux_residual(numComp); - bool converged_Well = true; - - // Finish computation - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) - { - 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 - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); - - if (std::isnan(well_flux_residual[phaseIdx])) { - OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); - } - - if (well_flux_residual[phaseIdx] > maxResidualAllowed) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); - } - } - - if ( terminal_output_ ) - { - // Only rank 0 does print to std::cout - if (iteration == 0) { - std::string msg; - msg = "Iter"; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); - msg += " W-FLUX(" + phaseName + ")"; + if (nan_residual_found) { + for (const auto& well : report.nan_residual_wells) { + OpmLog::debug("NaN residual found with phase " + well.phase_name + " for well " + well.well_name); } - OpmLog::note(msg); + OPM_THROW(Opm::NumericalProblem, "NaN residual found!"); } - - std::ostringstream ss; - const std::streamsize oprec = ss.precision(3); - const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); - ss << std::setw(4) << iteration; - 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()); } - return converged_Well; + + // checking too large residuals + { + bool too_large_residual_found = report.too_large_residual_found; + const auto& grid = ebosSimulator.gridManager().grid(); + int value = too_large_residual_found ? 1 : 0; + + too_large_residual_found = grid.comm().max(value); + if (too_large_residual_found) { + for (const auto& well : report.too_large_residual_wells) { + OpmLog::debug("Too large residual found with phase " + well.phase_name + " fow well " + well.well_name); + } + OPM_THROW(Opm::NumericalProblem, "Too large residual found!"); + } + } + + // checking convergence + bool converged_well = report.converged; + { + const auto& grid = ebosSimulator.gridManager().grid(); + int value = converged_well ? 1 : 0; + + converged_well = grid.comm().min(value); + } + + return converged_well; } @@ -1180,439 +536,13 @@ namespace Opm { void StandardWellsDense:: computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw) + const WellState& xw) const { if( ! localWellsActive() ) return ; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. - std::vector b_perf; - std::vector rsmax_perf; - std::vector rvmax_perf; - std::vector surf_dens_perf; - computePropertiesForWellConnectionPressures(ebosSimulator, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, cell_depths_, gravity_); - } - - - - - - template - void - StandardWellsDense:: - 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 - { - 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_; - 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]) { - rsmax_perf.resize(nperf); - rvmax_perf.resize(nperf); - } - - // Compute the average pressure in each well block - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - - const int cell_idx = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - - const double p_above = perf == wells().well_connpos[w] ? xw.bhp()[w] : xw.perfPress()[perf - 1]; - const double p_avg = (xw.perfPress()[perf] + p_above)/2; - const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); - - if (pu.phase_used[BlackoilPhases::Aqua]) { - 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 * numComp; - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; - - if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; - 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) { - const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); - double rv = 0.0; - if (gasrate > 0) { - rv = oilrate / gasrate; - } - rv = std::min(rv, rvmax_perf[perf]); - - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); - } - else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - - } else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; - 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 double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); - if (gasrate > 0) { - const double oilrate = std::abs(xw.wellRates()[oilpos_well]); - double rs = 0.0; - if (oilrate > 0) { - rs = gasrate / oilrate; - } - rs = std::min(rs, rsmax_perf[perf]); - b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - // Surface density. - for (int p = 0; p < pu.num_phases; ++p) { - surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); - } - - // We use cell values for solvent injector - if (has_solvent_) { - b_perf[numComp*perf + solventSaturationIdx] = intQuants.solventInverseFormationVolumeFactor().value(); - surf_dens_perf[numComp*perf + solventSaturationIdx] = intQuants.solventRefDensity(); - } - } - } - } - - - - - - template - void - StandardWellsDense:: - updateWellState(const BVector& dwells, - WellState& well_state) const - { - if( !localWellsActive() ) return; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - double dFLimit = dWellFractionMax(); - double dBHPLimit = dbhpMaxRel(); - std::vector xvar_well_old = well_state.wellSolutions(); - - for (int w = 0; w < nw; ++w) { - - // update the second and third well variable (The flux fractions) - std::vector F(np,0.0); - if (active_[ Water ]) { - const int sign2 = dwells[w][WFrac] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[w][WFrac]),dFLimit); - well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; - } - - if (active_[ Gas ]) { - const int sign3 = dwells[w][GFrac] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[w][GFrac]),dFLimit); - well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; - } - - if (has_solvent_) { - const int sign4 = dwells[w][SFrac] > 0 ? 1: -1; - const double dx4_limited = sign4 * std::min(std::abs(dwells[w][SFrac]),dFLimit); - well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited; - } - - assert(active_[ Oil ]); - F[Oil] = 1.0; - if (active_[ Water ]) { - F[Water] = well_state.wellSolutions()[WFrac*nw + w]; - F[Oil] -= F[Water]; - } - - if (active_[ Gas ]) { - F[Gas] = well_state.wellSolutions()[GFrac*nw + w]; - F[Oil] -= F[Gas]; - } - - double F_solvent = 0.0; - if (has_solvent_) { - F_solvent = well_state.wellSolutions()[SFrac*nw + w]; - F[Oil] -= F_solvent; - } - - if (active_[ Water ]) { - if (F[Water] < 0.0) { - if (active_[ Gas ]) { - F[Gas] /= (1.0 - F[Water]); - } - if (has_solvent_) { - F_solvent /= (1.0 - F[Water]); - } - F[Oil] /= (1.0 - F[Water]); - F[Water] = 0.0; - } - } - if (active_[ Gas ]) { - if (F[Gas] < 0.0) { - if (active_[ Water ]) { - F[Water] /= (1.0 - F[Gas]); - } - if (has_solvent_) { - F_solvent /= (1.0 - F[Gas]); - } - F[Oil] /= (1.0 - F[Gas]); - F[Gas] = 0.0; - } - } - if (F[Oil] < 0.0) { - if (active_[ Water ]) { - F[Water] /= (1.0 - F[Oil]); - } - if (active_[ Gas ]) { - F[Gas] /= (1.0 - F[Oil]); - } - if (has_solvent_) { - F_solvent /= (1.0 - F[Oil]); - } - F[Oil] = 0.0; - } - - if (active_[ Water ]) { - well_state.wellSolutions()[WFrac*nw + w] = F[Water]; - } - if (active_[ Gas ]) { - well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; - } - if(has_solvent_) { - well_state.wellSolutions()[SFrac*nw + w] = F_solvent; - } - - // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. - // More testing is needed to make sure this is correct for well groups and THP. - if (has_solvent_){ - F[Gas] += F_solvent; - } - - // The interpretation of the first well variable depends on the well control - const WellControls* wc = wells().ctrls[w]; - - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - const int current = well_state.currentControls()[w]; - const double target_rate = well_controls_iget_target(wc, current); - - std::vector g = {1,1,0.01}; - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - const double* distr = well_controls_iget_distr(wc, current); - for (int p = 0; p < np; ++p) { - if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value - F[p] /= distr[p]; - } else { - F[p] = 0.; - } - } - } else { - for (int p = 0; p < np; ++p) { - F[p] /= g[p]; - } - } - - switch (well_controls_iget_type(wc, current)) { - case THP: // The BHP and THP both uses the total rate as first well variable. - case BHP: - { - well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][XvarWell]; - - switch (wells().type[w]) { - case INJECTOR: - for (int p = 0; p < np; ++p) { - const double comp_frac = wells().comp_frac[np*w + p]; - well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[nw*XvarWell + w]; - } - break; - case PRODUCER: - for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = well_state.wellSolutions()[nw*XvarWell + w] * F[p]; - } - break; - } - - if (well_controls_iget_type(wc, current) == THP) { - - // Calculate bhp from thp control and well rates - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - // pick the density in the top layer - const int perf = wells().well_connpos[w]; - const double rho = well_perforation_densities_[perf]; - - if (well_type == INJECTOR) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - well_state.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - well_state.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - } - } - break; - case SURFACE_RATE: // Both rate controls use bhp as first well variable - case RESERVOIR_RATE: - { - const int sign1 = dwells[w][XvarWell] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[w][XvarWell]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit); - well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5); - well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w]; - - if (well_controls_iget_type(wc, current) == SURFACE_RATE) { - if (wells().type[w]==PRODUCER) { - - const double* distr = well_controls_iget_distr(wc, current); - - double F_target = 0.0; - for (int p = 0; p < np; ++p) { - F_target += distr[p] * F[p]; - } - for (int p = 0; p < np; ++p) { - well_state.wellRates()[np*w + p] = F[p] * target_rate / F_target; - } - } else { - - for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate; - } - } - } else { // RESERVOIR_RATE - for (int p = 0; p < np; ++p) { - well_state.wellRates()[np*w + p] = F[p] * target_rate; - } - } - } - break; - } // end of switch (well_controls_iget_type(wc, current)) - } // end of for (int w = 0; w < nw; ++w) - - - // for the wells having a THP constaint, we should update their thp value - // If it is under THP control, it will be set to be the target value. Otherwise, - // the thp value will be calculated based on the bhp value, assuming the bhp value is correctly calculated. - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - const int nwc = well_controls_get_num(wc); - // Looping over all controls until we find a THP constraint - int ctrl_index = 0; - for ( ; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(wc, ctrl_index) == THP) { - // the current control - const int current = well_state.currentControls()[w]; - // If under THP control at the moment - if (current == ctrl_index) { - const double thp_target = well_controls_iget_target(wc, current); - well_state.thp()[w] = thp_target; - } else { // otherwise we calculate the thp from the bhp value - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - const double alq = well_controls_iget_alq(wc, ctrl_index); - const int table_id = well_controls_iget_vfp(wc, ctrl_index); - - const WellType& well_type = wells().type[w]; - const int perf = wells().well_connpos[w]; //first perforation. - if (well_type == INJECTOR) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); - - const double bhp = well_state.bhp()[w]; - well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); - } else if (well_type == PRODUCER) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); - - const double bhp = well_state.bhp()[w]; - well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); - } else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - } - - // the THP control is found, we leave the loop now - break; - } - } // end of for loop for seaching THP constraints - - // no THP constraint found - if (ctrl_index == nwc) { // not finding a THP contstraints - well_state.thp()[w] = 0.0; - } - } // end of for (int w = 0; w < nw; ++w) + for (auto& well : well_container_) { + well->computeWellConnectionPressures(ebosSimulator, xw); + } } @@ -1630,97 +560,13 @@ namespace Opm { // we simply return. if( !wellsActive() ) return ; - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - // keeping a copy of the current controls, to see whether control changes later. - std::vector old_control_index(nw, 0); - for (int w = 0; w < nw; ++w) { - old_control_index[w] = xw.currentControls()[w]; - } - - // Find, for each well, if any constraints are broken. If so, - // switch control to first broken constraint. -#pragma omp parallel for schedule(dynamic) - for (int w = 0; w < nw; ++w) { - WellControls* wc = wells().ctrls[w]; - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - int current = xw.currentControls()[w]; - // Loop over all controls except the current one, and also - // skip any RESERVOIR_RATE controls, since we cannot - // handle those. - const int nwc = well_controls_get_num(wc); - int ctrl_index = 0; - for (; ctrl_index < nwc; ++ctrl_index) { - if (ctrl_index == current) { - // This is the currently used control, so it is - // used as an equation. So this is not used as an - // inequality constraint, and therefore skipped. - continue; - } - if (wellhelpers::constraintBroken( - xw.bhp(), xw.thp(), xw.wellRates(), - w, np, wells().type[w], wc, ctrl_index)) { - // ctrl_index will be the index of the broken constraint after the loop. - break; - } - } - if (ctrl_index != nwc) { - // Constraint number ctrl_index was broken, switch to it. - xw.currentControls()[w] = ctrl_index; - current = xw.currentControls()[w]; - well_controls_set_current( wc, current); - } - - // update whether well is under group control - if (wellCollection()->groupControlActive()) { - // get well node in the well collection - WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w])); - - // update whehter the well is under group control or individual control - if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) { - // under group control - well_node.setIndividualControl(false); - } else { - // individual control - well_node.setIndividualControl(true); - } - } - } - - // the new well control indices after all the related updates, - std::vector updated_control_index(nw, 0); - for (int w = 0; w < nw; ++w) { - updated_control_index[w] = xw.currentControls()[w]; - } - - // checking whether control changed wellhelpers::WellSwitchingLogger logger; - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - if (updated_control_index[w] != old_control_index[w]) { - logger.wellSwitched(wells().name[w], - well_controls_iget_type(wc, old_control_index[w]), - well_controls_iget_type(wc, updated_control_index[w])); - } - if (updated_control_index[w] != old_control_index[w] || well_collection_->groupControlActive()) { - updateWellStateWithTarget(wc, updated_control_index[w], w, xw); - } + for (const auto& well : well_container_) { + well->updateWellControl(xw, logger); } - // upate the well targets following group controls - // it will not change the control mode, only update the targets - if (wellCollection()->groupControlActive()) { - applyVREPGroupControl(xw); - wellCollection()->updateWellTargets(xw.wellRates()); - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - updateWellStateWithTarget(wc, updated_control_index[w], w, xw); - } - } + updateGroupControls(xw); } @@ -1736,138 +582,9 @@ namespace Opm { const WellState& well_state, DynamicListEconLimited& list_econ_limited) const { - // With no wells (on process) wells_struct is a null pointer - const int nw = (wells_struct)? wells_struct->number_of_wells : 0; - - for (int w = 0; w < nw; ++w) { - // flag to check if the mim oil/gas rate limit is violated - bool rate_limit_violated = false; - const std::string& well_name = wells_struct->name[w]; - const Well* well_ecl = schedule.getWell(well_name); - const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step); - - // economic limits only apply for production wells. - if (wells_struct->type[w] != PRODUCER) { - continue; - } - - // if no limit is effective here, then continue to the next well - if ( !econ_production_limits.onAnyEffectiveLimit() ) { - continue; - } - // for the moment, we only handle rate limits, not handling potential limits - // the potential limits should not be difficult to add - const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); - if (quantity_limit == WellEcon::POTN) { - const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") - + std::string("All the limits will be evaluated based on RATE. "); - OpmLog::warning("NOT_SUPPORTING_POTN", msg); - } - - const WellMapType& well_map = well_state.wellMap(); - const typename WellMapType::const_iterator i_well = well_map.find(well_name); - assert(i_well != well_map.end()); // should always be found? - const WellMapEntryType& map_entry = i_well->second; - const int well_number = map_entry[0]; - - if (econ_production_limits.onAnyRateLimit()) { - rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number); - } - - if (rate_limit_violated) { - if (econ_production_limits.endRun()) { - const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") - + std::string("the program will keep running after ") + well_name + std::string(" is closed"); - OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); - } - - if (econ_production_limits.validFollowonWell()) { - OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); - } - - if (well_ecl->getAutomaticShutIn()) { - list_econ_limited.addShutWell(well_name); - const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); - OpmLog::info(msg); - } else { - list_econ_limited.addStoppedWell(well_name); - const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); - OpmLog::info(msg); - } - // the well is closed, not need to check other limits - continue; - } - - // checking for ratio related limits, mostly all kinds of ratio. - bool ratio_limits_violated = false; - RatioCheckTuple ratio_check_return; - - if (econ_production_limits.onAnyRatioLimit()) { - ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry); - ratio_limits_violated = std::get<0>(ratio_check_return); - } - - if (ratio_limits_violated) { - const bool last_connection = std::get<1>(ratio_check_return); - const int worst_offending_connection = std::get<2>(ratio_check_return); - - const int perf_start = map_entry[1]; - - assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2])); - - const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection]; - list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); - const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") - + well_name + std::string(" will be closed due to economic limit"); - OpmLog::info(msg); - - if (last_connection) { - list_econ_limited.addShutWell(well_name); - const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); - OpmLog::info(msg2); - } - } - - } // for (int w = 0; w < nw; ++w) - } - - - - - - template - void - StandardWellsDense:: - computeWellConnectionDensitesPressures(const WellState& xw, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf, - const std::vector& depth_perf, - const double grav) - { - // Compute densities - const int nperf = depth_perf.size(); - const int numComponent = b_perf.size() / nperf; - const int np = wells().number_of_phases; - std::vector perfRates(b_perf.size(),0.0); - for (int perf = 0; perf < nperf; ++perf) { - for (int phase = 0; phase < np; ++phase) { - perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase]; - } - if(has_solvent_) { - perfRates[perf*numComponent + solventSaturationIdx] = xw.perfRateSolvent()[perf]; - } + for (const auto& well : well_container_) { + well->updateListEconLimited(well_state, list_econ_limited); } - well_perforation_densities_ = - WellDensitySegmented::computeConnectionDensities( - wells(), phase_usage_, perfRates, - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // Compute pressure deltas - well_perforation_pressure_diffs_ = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), depth_perf, well_perforation_densities_, grav); } @@ -1881,50 +598,21 @@ namespace Opm { const WellState& well_state, std::vector& well_potentials) const { + updatePrimaryVariables(well_state); + computeWellConnectionPressures(ebosSimulator, well_state); + + // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials + // TODO: for computeWellPotentials, no derivative is required actually + initPrimaryVariablesEvaluation(); // number of wells and phases - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - + const int nw = number_of_wells_; + const int np = number_of_phases_; well_potentials.resize(nw * np, 0.0); for (int w = 0; w < nw; ++w) { - - // get the bhp value based on the bhp constraints - const double bhp = mostStrictBhpFromBhpLimits(w); - - // does the well have a THP related constraint? - const bool has_thp_control = wellHasTHPConstraints(w); - - std::vector potentials(np); - - if ( !has_thp_control ) { - - assert(std::abs(bhp) != std::numeric_limits::max()); - - computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials); - - } else { // the well has a THP related constraint - // checking whether a well is newly added, it only happens at the beginning of the report step - if ( !well_state.isNewWell(w) ) { - for (int p = 0; p < np; ++p) { - // This is dangerous for new added well - // since we are not handling the initialization correctly for now - potentials[p] = well_state.wellRates()[w * np + p]; - } - } else { - // We need to generate a reasonable rates to start the iteration process - computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials); - for (double& value : potentials) { - // make the value a little safer in case the BHP limits are default ones - // TODO: a better way should be a better rescaling based on the investigation of the VFP table. - const double rate_safety_scaling_factor = 0.00001; - value *= rate_safety_scaling_factor; - } - } - - potentials = computeWellPotentialWithTHP(ebosSimulator, w, bhp, potentials); - } + std::vector potentials; + well_container_[w]->computeWellPotentials(ebosSimulator, well_state, potentials); // putting the sucessfully calculated potentials to the well_potentials for (int p = 0; p < np; ++p) { @@ -1943,7 +631,7 @@ namespace Opm { prepareTimeStep(const Simulator& ebos_simulator, WellState& well_state) { - const int nw = wells().number_of_wells; + const int nw = number_of_wells_; for (int w = 0; w < nw; ++w) { // after restarting, the well_controls can be modified while // the well_state still uses the old control index @@ -1951,8 +639,8 @@ namespace Opm { resetWellControlFromState(well_state); if (wellCollection()->groupControlActive()) { - WellControls* wc = wells().ctrls[w]; - WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w])); + WellControls* wc = well_container_[w]->wellControls(); + WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name()); // handling the situation that wells do not have a valid control // it happens the well specified with GRUP and restarting due to non-convergencing @@ -1983,15 +671,13 @@ namespace Opm { if (well_collection_->requireWellPotentials()) { // calculate the well potentials - setWellVariables(well_state); - computeWellConnectionPressures(ebos_simulator, well_state); - - // To store well potentials for each well std::vector well_potentials; computeWellPotentials(ebos_simulator, well_state, well_potentials); // update/setup guide rates for each well based on the well_potentials - well_collection_->setGuideRatesWithPotentials(wellsPointer(), phase_usage_, well_potentials); + // TODO: this is one of two places that still need Wells struct. In this function, only the well names + // well types are used, probably the order of the wells to locate the correct values in well_potentials. + well_collection_->setGuideRatesWithPotentials(wells_, phase_usage_, well_potentials); } applyVREPGroupControl(well_state); @@ -2005,10 +691,12 @@ namespace Opm { // since the controls are all updated, we should update well_state accordingly for (int w = 0; w < nw; ++w) { - WellControls* wc = wells().ctrls[w]; + WellControls* wc = well_container_[w]->wellControls(); const int control = well_controls_get_current(wc); well_state.currentControls()[w] = control; - updateWellStateWithTarget(wc, control, w, well_state); + // TODO: for VFP control, the perf_densities are still zero here, investigate better + // way to handle it later. + well_container_[w]->updateWellStateWithTarget(control, well_state); // The wells are not considered to be newly added // for next time step @@ -2033,17 +721,6 @@ namespace Opm { - template - const std::vector& - StandardWellsDense:: - wellPerfEfficiencyFactors() const - { - return well_perforation_efficiency_factors_; - } - - - - template void @@ -2054,18 +731,15 @@ namespace Opm { return; } - const int nw = wells().number_of_wells; + const int nw = number_of_wells_; for (int w = 0; w < nw; ++w) { - const std::string well_name = wells().name[w]; + const std::string well_name = well_container_[w]->name(); const WellNode& well_node = wellCollection()->findWellNode(well_name); const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor(); - // assign the efficiency factor to each perforation related. - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) { - well_perforation_efficiency_factors_[perf] = well_efficiency_factor; - } + well_container_[w]->setWellEfficiencyFactor(well_efficiency_factor); } } @@ -2088,8 +762,8 @@ namespace Opm { // We only store the conversion coefficients for all the injection wells. // Later, more delicate model will be implemented here. // And for the moment, group control can only work for serial running. - const int nw = well_state.numWells(); - const int np = well_state.numPhases(); + const int nw = numWells(); + const int np = numPhases(); // we calculate the voidage rate for each well, that means the sum of all the phases. well_voidage_rates.resize(nw, 0); @@ -2100,7 +774,7 @@ namespace Opm { std::vector convert_coeff(np, 1.0); for (int w = 0; w < nw; ++w) { - const bool is_producer = wells().type[w] == PRODUCER; + const bool is_producer = well_container_[w]->wellType() == PRODUCER; // not sure necessary to change all the value to be positive if (is_producer) { @@ -2150,7 +824,7 @@ namespace Opm { const int well_index = well_node->selfIndex(); well_state.currentControls()[well_index] = well_node->groupControlIndex(); - WellControls* wc = wells().ctrls[well_index]; + WellControls* wc = well_container_[well_index]->wellControls(); well_controls_set_current(wc, well_node->groupControlIndex()); } } @@ -2161,740 +835,41 @@ namespace Opm { - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - getBhp(const int wellIdx) const { - const WellControls* wc = wells().ctrls[wellIdx]; - if (well_controls_get_current_type(wc) == BHP) { - EvalWell bhp = 0.0; - const double target_rate = well_controls_get_current_target(wc); - bhp.setValue(target_rate); - return bhp; - } else if (well_controls_get_current_type(wc) == THP) { - const int control = well_controls_get_current(wc); - const double thp = well_controls_get_current_target(wc); - const double alq = well_controls_iget_alq(wc, control); - const int table_id = well_controls_iget_vfp(wc, control); - EvalWell aqua = 0.0; - EvalWell liquid = 0.0; - EvalWell vapour = 0.0; - EvalWell bhp = 0.0; - double vfp_ref_depth = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = getQs(wellIdx, pu.phase_pos[ Water]); - } - if (active_[ Oil ]) { - liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); - } - if (active_[ Gas ]) { - vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); - } - if (wells().type[wellIdx] == INJECTOR) { - bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); - vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); - } else { - bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); - vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); - } - - // pick the density in the top layer - const int perf = wells().well_connpos[wellIdx]; - const double rho = well_perforation_densities_[perf]; - const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); - bhp -= dp; - return bhp; - } - - const int nw = wells().number_of_wells; - return wellVariables_[nw*XvarWell + wellIdx]; - } - - - - - - 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 < numComponents()); - const int nw = wells().number_of_wells; - const auto pu = phase_usage_; - 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) { - if (has_solvent_ ) { - double comp_frac = 0.0; - if (has_solvent_ && compIdx == solventSaturationIdx) { // solvent - comp_frac = wells().comp_frac[np*wellIdx + pu.phase_pos[ Gas ]] * wsolvent(wellIdx); - } else if (compIdx == pu.phase_pos[ Gas ]) { - comp_frac = wells().comp_frac[np*wellIdx + compIdx] * (1.0 - wsolvent(wellIdx)); - } else { - comp_frac = wells().comp_frac[np*wellIdx + compIdx]; - } - if (comp_frac == 0.0) { - return qs; //zero - } - - if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { - return comp_frac * wellVariables_[nw*XvarWell + wellIdx]; - } - - qs.setValue(comp_frac * target_rate); - return qs; - } - const double comp_frac = wells().comp_frac[np*wellIdx + compIdx]; - if (comp_frac == 0.0) { - return qs; - } - - if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { - return wellVariables_[nw*XvarWell + wellIdx]; - } - qs.setValue(target_rate); - return qs; - } - - // Producers - if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { - return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx, compIdx); - } - - if (well_controls_get_current_type(wc) == SURFACE_RATE) { - // checking how many phases are included in the rate control - // to decide wheter it is a single phase rate control or not - const double* distr = well_controls_get_current_distr(wc); - int num_phases_under_rate_control = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - num_phases_under_rate_control += 1; - } - } - - // there should be at least one phase involved - assert(num_phases_under_rate_control > 0); - - // when it is a single phase rate limit - if (num_phases_under_rate_control == 1) { - - // looking for the phase under control - int phase_under_control = -1; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - phase_under_control = phase; - break; - } - } - - assert(phase_under_control >= 0); - - EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(wellIdx, phase_under_control); - if (has_solvent_ && phase_under_control == Gas) { - // for GRAT controlled wells solvent is included in the target - wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventSaturationIdx); - } - - if (compIdx == phase_under_control) { - if (has_solvent_ && phase_under_control == Gas) { - qs.setValue(target_rate * wellVolumeFractionScaled(wellIdx, Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); - return qs; - } - qs.setValue(target_rate); - return qs; - } - - // TODO: not sure why the single phase under control will have near zero fraction - const double eps = 1e-6; - if (wellVolumeFractionScaledPhaseUnderControl < eps) { - return qs; - } - return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaledPhaseUnderControl); - } - - // when it is a combined two phase rate limit, such like LRAT - // we neec to calculate the rate for the certain phase - if (num_phases_under_rate_control == 2) { - EvalWell combined_volume_fraction = 0.; - for (int p = 0; p < np; ++p) { - if (distr[p] == 1.0) { - combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p); - } - } - 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, compIdx); - } - } else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - // ReservoirRate - return target_rate * wellVolumeFractionScaled(wellIdx, compIdx); - } else { - OPM_THROW(std::logic_error, "Unknown control type for well " << wells().name[wellIdx]); - } - - // avoid warning of condition reaches end of non-void function - return qs; - } - - - - - - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - wellVolumeFraction(const int wellIdx, const int compIdx) const - { - const int nw = wells().number_of_wells; - if (compIdx == Water) { - return wellVariables_[WFrac * nw + wellIdx]; - } - - if (compIdx == Gas) { - return wellVariables_[GFrac * nw + wellIdx]; - } - - if (has_solvent_ && compIdx == solventSaturationIdx) { - return wellVariables_[SFrac * nw + wellIdx]; - } - - // Oil fraction - EvalWell well_fraction = 1.0; - if (active_[Water]) { - well_fraction -= wellVariables_[WFrac * nw + wellIdx]; - } - - if (active_[Gas]) { - well_fraction -= wellVariables_[GFrac * nw + wellIdx]; - } - if (has_solvent_) { - well_fraction -= wellVariables_[SFrac * nw + wellIdx]; - } - return well_fraction; - } - - - - - - 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) { - - if (has_solvent_ && compIdx == solventSaturationIdx) { - return wellVolumeFraction(wellIdx, compIdx); - } - const double* distr = well_controls_get_current_distr(wc); - 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, compIdx); - } - } - std::vector g = {1,1,0.01,0.01}; - return (wellVolumeFraction(wellIdx, compIdx) / g[compIdx]); - } - - - - - - template - typename StandardWellsDense::EvalWell - StandardWellsDense:: - wellSurfaceVolumeFraction(const int well_index, const int compIdx) const - { - EvalWell sum_volume_fraction_scaled = 0.; - 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, compIdx) / sum_volume_fraction_scaled; - } - - - - - - template - bool - 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(); - - if (econ_production_limits.onMinOilRate()) { - assert(active_[Oil]); - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double min_oil_rate = econ_production_limits.minOilRate(); - if (std::abs(oil_rate) < min_oil_rate) { - return true; - } - } - - if (econ_production_limits.onMinGasRate() ) { - assert(active_[Gas]); - const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; - const double min_gas_rate = econ_production_limits.minGasRate(); - if (std::abs(gas_rate) < min_gas_rate) { - return true; - } - } - - if (econ_production_limits.onMinLiquidRate() ) { - assert(active_[Oil]); - assert(active_[Water]); - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; - const double liquid_rate = oil_rate + water_rate; - const double min_liquid_rate = econ_production_limits.minLiquidRate(); - if (std::abs(liquid_rate) < min_liquid_rate) { - return true; - } - } - - if (econ_production_limits.onMinReservoirFluidRate()) { - OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); - } - - return false; - } - - - - - - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: - checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const - { - // TODO: not sure how to define the worst-offending connection when more than one - // ratio related limit is violated. - // The defintion used here is that we define the violation extent based on the - // ratio between the value and the corresponding limit. - // For each violated limit, we decide the worst-offending connection separately. - // Among the worst-offending connections, we use the one has the biggest violation - // extent. - - bool any_limit_violated = false; - bool last_connection = false; - int worst_offending_connection = INVALIDCONNECTION; - double violation_extent = -1.0; - - if (econ_production_limits.onMaxWaterCut()) { - const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry); - bool water_cut_violated = std::get<0>(water_cut_return); - if (water_cut_violated) { - any_limit_violated = true; - const double violation_extent_water_cut = std::get<3>(water_cut_return); - if (violation_extent_water_cut > violation_extent) { - violation_extent = violation_extent_water_cut; - worst_offending_connection = std::get<2>(water_cut_return); - last_connection = std::get<1>(water_cut_return); - } - } - } - - if (econ_production_limits.onMaxGasOilRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); - } - - if (econ_production_limits.onMaxWaterGasRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); - } - - if (econ_production_limits.onMaxGasLiquidRatio()) { - OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); - } - - if (any_limit_violated) { - assert(worst_offending_connection >=0); - assert(violation_extent > 1.); - } - - return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - - - - - template - typename StandardWellsDense::RatioCheckTuple - StandardWellsDense:: - checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, - const WellState& well_state, - const WellMapEntryType& map_entry) const - { - bool water_cut_limit_violated = false; - int worst_offending_connection = INVALIDCONNECTION; - bool last_connection = false; - double violation_extent = -1.0; - - const int np = well_state.numPhases(); - const Opm::PhaseUsage& pu = phase_usage_; - const int well_number = map_entry[0]; - - assert(active_[Oil]); - assert(active_[Water]); - - const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; - const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; - const double liquid_rate = oil_rate + water_rate; - double water_cut; - if (std::abs(liquid_rate) != 0.) { - water_cut = water_rate / liquid_rate; - } else { - water_cut = 0.0; - } - - const double max_water_cut_limit = econ_production_limits.maxWaterCut(); - if (water_cut > max_water_cut_limit) { - water_cut_limit_violated = true; - } - - if (water_cut_limit_violated) { - // need to handle the worst_offending_connection - const int perf_start = map_entry[1]; - const int perf_number = map_entry[2]; - - std::vector water_cut_perf(perf_number); - for (int perf = 0; perf < perf_number; ++perf) { - const int i_perf = perf_start + perf; - const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; - const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; - const double liquid_perf_rate = oil_perf_rate + water_perf_rate; - if (std::abs(liquid_perf_rate) != 0.) { - water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; - } else { - water_cut_perf[perf] = 0.; - } - } - - last_connection = (perf_number == 1); - if (last_connection) { - worst_offending_connection = 0; - violation_extent = water_cut_perf[0] / max_water_cut_limit; - return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - double max_water_cut_perf = 0.; - for (int perf = 0; perf < perf_number; ++perf) { - if (water_cut_perf[perf] > max_water_cut_perf) { - worst_offending_connection = perf; - max_water_cut_perf = water_cut_perf[perf]; - } - } - - assert(max_water_cut_perf != 0.); - assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); - - violation_extent = max_water_cut_perf / max_water_cut_limit; - } - - return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); - } - - - - - template void StandardWellsDense:: - updateWellStateWithTarget(const WellControls* wc, - const int current, - const int well_index, - WellState& xw) const + updateGroupControls(WellState& well_state) const { - // number of phases - const int np = wells().number_of_phases; - // Updating well state and primary variables. - // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE - const double target = well_controls_iget_target(wc, current); - const double* distr = well_controls_iget_distr(wc, current); - switch (well_controls_iget_type(wc, current)) { - case BHP: - xw.bhp()[well_index] = target; - // TODO: similar to the way below to handle THP - // we should not something related to thp here when there is thp constraint - break; - case THP: { - xw.thp()[well_index] = target; + if (wellCollection()->groupControlActive()) { + for (int w = 0; w < number_of_wells_; ++w) { + // update whether well is under group control + // get well node in the well collection + WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name()); - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[well_index]; - - // pick the density in the top layer - const int perf = wells().well_connpos[well_index]; - const double rho = well_perforation_densities_[perf]; - - if (well_type == INJECTOR) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - rho, gravity_); - - xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - break; - } - - case RESERVOIR_RATE: // intentional fall-through - case SURFACE_RATE: - // checking the number of the phases under control - int numPhasesWithTargetsUnderThisControl = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - numPhasesWithTargetsUnderThisControl += 1; + // update whehter the well is under group control or individual control + const int current = well_state.currentControls()[w]; + if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) { + // under group control + well_node.setIndividualControl(false); + } else { + // individual control + well_node.setIndividualControl(true); } } - assert(numPhasesWithTargetsUnderThisControl > 0); + applyVREPGroupControl(well_state); + // upate the well targets following group controls + // it will not change the control mode, only update the targets + wellCollection()->updateWellTargets(well_state.wellRates()); - const WellType& well_type = wells().type[well_index]; - if (well_type == INJECTOR) { - // assign target value as initial guess for injectors - // only handles single phase control at the moment - assert(numPhasesWithTargetsUnderThisControl == 1); - - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.) { - xw.wellRates()[np*well_index + phase] = target / distr[phase]; - } else { - xw.wellRates()[np * well_index + phase] = 0.; - } - } - } else if (well_type == PRODUCER) { - // update the rates of phases under control based on the target, - // and also update rates of phases not under control to keep the rate ratio, - // assuming the mobility ratio does not change for the production wells - double original_rates_under_phase_control = 0.0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase]; - } - } - - if (original_rates_under_phase_control != 0.0 ) { - double scaling_factor = target / original_rates_under_phase_control; - - for (int phase = 0; phase < np; ++phase) { - xw.wellRates()[np * well_index + phase] *= scaling_factor; - } - } else { // scaling factor is not well defied when original_rates_under_phase_control is zero - // separating targets equally between phases under control - const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase]; - } else { - // this only happens for SURFACE_RATE control - xw.wellRates()[np * well_index + phase] = target_rate_divided; - } - } - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - break; - } // end of switch - - - std::vector g = {1.0, 1.0, 0.01}; - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - for (int phase = 0; phase < np; ++phase) { - g[phase] = distr[phase]; - } - } - - // the number of wells - const int nw = wells().number_of_wells; - - switch (well_controls_iget_type(wc, current)) { - case THP: - case BHP: { - const WellType& well_type = wells().type[well_index]; - xw.wellSolutions()[nw*XvarWell + well_index] = 0.0; - if (well_type == INJECTOR) { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * wells().comp_frac[np*well_index + p]; - } - } else { - for (int p = 0; p < np; ++p) { - xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p]; - } - } - break; - } - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index]; - break; - } // end of switch - - double tot_well_rate = 0.0; - for (int p = 0; p < np; ++p) { - tot_well_rate += g[p] * xw.wellRates()[np*well_index + p]; - } - if(std::abs(tot_well_rate) > 0) { - if (active_[ Water ]) { - xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate; - } - if (active_[ Gas ]) { - xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (xw.wellRates()[np*well_index + Gas] - xw.solventWellRate(well_index)) / tot_well_rate ; - } - if (has_solvent_) { - xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * xw.solventWellRate(well_index) / tot_well_rate ; - } - } else { - const WellType& well_type = wells().type[well_index]; - if (well_type == INJECTOR) { - // only single phase injection handled - if (active_[Water]) { - if (distr[Water] > 0.0) { - xw.wellSolutions()[WFrac * nw + well_index] = 1.0; - } else { - xw.wellSolutions()[WFrac * nw + well_index] = 0.0; - } - } - - if (active_[Gas]) { - if (distr[Gas] > 0.0) { - xw.wellSolutions()[GFrac * nw + well_index] = 1.0 - wsolvent(well_index); - if (has_solvent_) { - xw.wellSolutions()[SFrac * nw + well_index] = wsolvent(well_index); - } - } else { - xw.wellSolutions()[GFrac * nw + well_index] = 0.0; - } - } - - // TODO: it is possible to leave injector as a oil well, - // when F_w and F_g both equals to zero, not sure under what kind of circumstance - // this will happen. - } else if (well_type == PRODUCER) { // producers - if (active_[Water]) { - xw.wellSolutions()[WFrac * nw + well_index] = 1.0 / np; - } - if (active_[Gas]) { - xw.wellSolutions()[GFrac * nw + well_index] = 1.0 / np; - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - } - - } - - - - - - template - bool - StandardWellsDense:: - wellHasTHPConstraints(const int well_index) const - { - const WellControls* well_control = wells().ctrls[well_index]; - const int nwc = well_controls_get_num(well_control); - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(well_control, ctrl_index) == THP) { - return true; - } - } - return false; - } - - - - - - template - void - StandardWellsDense:: - computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, - const int well_index, - std::vector& well_flux) const - { - const int np = wells().number_of_phases; - const int numComp = numComponents(); - well_flux.resize(np, 0.0); - - const bool allow_cf = allow_cross_flow(well_index, ebosSimulator); - for (int perf = wells().well_connpos[well_index]; perf < wells().well_connpos[well_index + 1]; ++perf) { - const int cell_index = wells().well_cells[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); - // flux for each perforation - std::vector cq_s(numComp, 0.0); - std::vector mob(numComp, 0.0); - getMobility(ebosSimulator, well_index, perf, cell_index, mob); - computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp, - wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - - for(int p = 0; p < np; ++p) { - well_flux[p] += cq_s[p].value(); + for (int w = 0; w < number_of_wells_; ++w) { + // TODO: check whether we need current argument in updateWellStateWithTarget + // maybe there is some circumstances that the current is different from the one + // in the WellState. + // while probalby, the current argument can be removed + const int current = well_state.currentControls()[w]; + well_container_[w]->updateWellStateWithTarget(current, well_state); } } } @@ -2903,255 +878,6 @@ namespace Opm { - - template - double - StandardWellsDense:: - mostStrictBhpFromBhpLimits(const int well_index) const - { - double bhp; - - // type of the well, INJECTOR or PRODUCER - const WellType& well_type = wells().type[well_index]; - // initial bhp value, making the value not usable - switch(well_type) { - case INJECTOR: - bhp = std::numeric_limits::max(); - break; - case PRODUCER: - bhp = -std::numeric_limits::max(); - break; - default: - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[well_index]); - } - - // the well controls - const WellControls* well_control = wells().ctrls[well_index]; - // The number of the well controls/constraints - const int nwc = well_controls_get_num(well_control); - - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - // finding a BHP constraint - if (well_controls_iget_type(well_control, ctrl_index) == BHP) { - // get the bhp constraint value, it should always be postive assummingly - const double bhp_target = well_controls_iget_target(well_control, ctrl_index); - - switch(well_type) { - case INJECTOR: // using the lower bhp contraint from Injectors - if (bhp_target < bhp) { - bhp = bhp_target; - } - break; - case PRODUCER: - if (bhp_target > bhp) { - bhp = bhp_target; - } - break; - default: - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[well_index]); - } // end of switch - } - } - - return bhp; - } - - - - - - template - std::vector - StandardWellsDense:: - computeWellPotentialWithTHP(const Simulator& ebosSimulator, - const int well_index, - const double initial_bhp, // bhp from BHP constraints - const std::vector& initial_potential) const - { - // TODO: pay attention to the situation that finally the potential is calculated based on the bhp control - // TODO: should we consider the bhp constraints during the iterative process? - const int np = wells().number_of_phases; - - assert( np == int(initial_potential.size()) ); - - std::vector potentials = initial_potential; - std::vector old_potentials = potentials; // keeping track of the old potentials - - double bhp = initial_bhp; - double old_bhp = bhp; - - bool converged = false; - const int max_iteration = 1000; - const double bhp_tolerance = 1000.; // 1000 pascal - - int iteration = 0; - - while ( !converged && iteration < max_iteration ) { - // for each iteration, we calculate the bhp based on the rates/potentials with thp constraints - // with considering the bhp value from the bhp limits. At the beginning of each iteration, - // we initialize the bhp to be the bhp value from the bhp limits. Then based on the bhp values calculated - // from the thp constraints, we decide the effective bhp value for well potential calculation. - bhp = initial_bhp; - - // the well controls - const WellControls* well_control = wells().ctrls[well_index]; - // The number of the well controls/constraints - const int nwc = well_controls_get_num(well_control); - - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(well_control, ctrl_index) == THP) { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - const Opm::PhaseUsage& pu = phase_usage_; - - if (active_[ Water ]) { - aqua = potentials[pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = potentials[pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = potentials[pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(well_control, ctrl_index); - const double thp = well_controls_iget_target(well_control, ctrl_index); - const double alq = well_controls_iget_alq(well_control, ctrl_index); - - // Calculating the BHP value based on THP - // TODO: check whether it is always correct to do calculation based on the depth of the first perforation. - const int first_perf = wells().well_connpos[well_index]; //first perforation - - const WellType& well_type = wells().type[well_index]; - if (well_type == INJECTOR) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[first_perf], gravity_); - const double bhp_calculated = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - // apply the strictest of the bhp controlls i.e. smallest bhp for injectors - if (bhp_calculated < bhp) { - bhp = bhp_calculated; - } - } - else if (well_type == PRODUCER) { - const double dp = wellhelpers::computeHydrostaticCorrection( - wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[first_perf], gravity_); - const double bhp_calculated = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - // apply the strictest of the bhp controlls i.e. largest bhp for producers - if (bhp_calculated > bhp) { - bhp = bhp_calculated; - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - } - } - - // there should be always some available bhp/thp constraints there - if (std::isinf(bhp) || std::isnan(bhp)) { - OPM_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << wells().name[well_index]); - } - - converged = std::abs(old_bhp - bhp) < bhp_tolerance; - - computeWellRatesWithBhp(ebosSimulator, bhp, well_index, potentials); - - // checking whether the potentials have valid values - for (const double value : potentials) { - if (std::isinf(value) || std::isnan(value)) { - OPM_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << wells().name[well_index]); - } - } - - if (!converged) { - old_bhp = bhp; - for (int p = 0; p < np; ++p) { - // TODO: improve the interpolation, will it always be valid with the way below? - // TODO: finding better paramters, better iteration strategy for better convergence rate. - const double potential_update_damping_factor = 0.001; - potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p]; - old_potentials[p] = potentials[p]; - } - } - - ++iteration; - } - - if (!converged) { - OPM_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << wells().name[well_index]); - } - - return potentials; - } - - template - double - StandardWellsDense:: - wsolvent(const int well_index) const { - - if (!has_solvent_) { - return 0.0; - } - - // loop over all wells until we find the well with the matching name - for (const auto& well : wells_ecl_) { - if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) { - continue; - } - - WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_); - if (injection.injectorType == WellInjector::GAS) { - - double solventFraction = well->getSolventFraction(current_timeIdx_); - - // Look until we find the correct well - if (well->name() == wells().name[well_index]) { - return solventFraction; - } - } - } - // we didn't find it return 0; - assert(false); - return 0.0; - } - - - template - double - StandardWellsDense:: - wpolymer(const int well_index) const { - - if (!has_polymer_) { - return 0.0; - } - - // loop over all wells until we find the well with the matching name - for (const auto& well : wells_ecl_) { - if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) { - continue; - } - - WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_); - WellPolymerProperties polymer = well->getPolymerProperties(current_timeIdx_); - if (injection.injectorType == WellInjector::WATER) { - - double polymerFraction = polymer.m_polymerConcentration; - - // Look until we find the correct well - if (well->name() == wells().name[well_index]) { - return polymerFraction; - } - } - } - // we didn't find it return 0; - assert(false); - return 0.0; - } - - template void StandardWellsDense:: @@ -3175,115 +901,78 @@ namespace Opm { StandardWellsDense:: computeRepRadiusPerfLength(const Grid& grid) { - // TODO, the function does not work for parallel running // to be fixed later. - int number_of_cells = Opm::UgGridHelpers::numCells(grid); const int* global_cell = Opm::UgGridHelpers::globalCell(grid); - const int* cart_dims = Opm::UgGridHelpers::cartDims(grid); - auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid); - auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid); - - if (wells_ecl_.size() == 0) { - OPM_MESSAGE("No wells specified in Schedule section, " - "initializing no wells"); - return; - } - - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - - const size_t timeStep = current_timeIdx_; - - wells_rep_radius_.clear(); - wells_perf_length_.clear(); - wells_bore_diameter_.clear(); - - wells_rep_radius_.reserve(nperf); - wells_perf_length_.reserve(nperf); - wells_bore_diameter_.reserve(nperf); std::map cartesian_to_compressed; - - setupCompressedToCartesian(global_cell, number_of_cells, + setupCompressedToCartesian(global_cell, number_of_cells_, cartesian_to_compressed); - int well_index = 0; - - for (auto wellIter= wells_ecl_.begin(); wellIter != wells_ecl_.end(); ++wellIter) { - const auto* well = (*wellIter); - - if (well->getStatus(timeStep) == WellCommon::SHUT) { - continue; - } - { // COMPDAT handling - const auto& completionSet = well->getCompletions(timeStep); - for (size_t c=0; c::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); - if (cgit == cartesian_to_compressed.end()) { - OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' - << k << " not found in grid (well = " << well->name() << ')'); - } - int cell = cgit->second; - - { - double radius = 0.5*completion.getDiameter(); - if (radius <= 0.0) { - radius = 0.5*unit::feet; - OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); - } - - const std::array cubical = - WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell); - - WellCompletion::DirectionEnum direction = completion.getDirection(); - - double re; // area equivalent radius of the grid block - double perf_length; // the length of the well perforation - - switch (direction) { - case Opm::WellCompletion::DirectionEnum::X: - re = std::sqrt(cubical[1] * cubical[2] / M_PI); - perf_length = cubical[0]; - break; - case Opm::WellCompletion::DirectionEnum::Y: - re = std::sqrt(cubical[0] * cubical[2] / M_PI); - perf_length = cubical[1]; - break; - case Opm::WellCompletion::DirectionEnum::Z: - re = std::sqrt(cubical[0] * cubical[1] / M_PI); - perf_length = cubical[2]; - break; - default: - OPM_THROW(std::runtime_error, " Dirtecion of well is not supported "); - } - - double repR = std::sqrt(re * radius); - wells_rep_radius_.push_back(repR); - wells_perf_length_.push_back(perf_length); - wells_bore_diameter_.push_back(2. * radius); - } - } else { - if (completion.getState() != WellCompletion::SHUT) { - OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled"); - } - } - - } - } - well_index++; + for (const auto& well : well_container_) { + well->computeRepRadiusPerfLength(grid, cartesian_to_compressed); } } + + template + void + StandardWellsDense:: + computeAverageFormationFactor(Simulator& ebosSimulator, + std::vector& B_avg) const + { + const int np = numPhases(); + + const auto& grid = ebosSimulator.gridManager().grid(); + const auto& gridView = grid.leafGridView(); + ElementContext elemCtx(ebosSimulator); + const auto& elemEndIt = gridView.template end(); + + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; ++elemIt) + { + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) + { + auto& B = B_avg[ phaseIdx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); + + B += 1 / fs.invB(ebosPhaseIdx).value(); + } + if (has_solvent_) { + auto& B = B_avg[solventSaturationIdx]; + B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); + } + } + + // compute global average + grid.comm().sum(B_avg.data(), B_avg.size()); + for(auto& bval: B_avg) + { + bval/=global_nc_; + } + } + + + + + + template + void + StandardWellsDense:: + updatePrimaryVariables(const WellState& well_state) const + { + for (const auto& well : well_container_) { + well->updatePrimaryVariables(well_state); + } + } + } // namespace Opm diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp index c26168ef3..8bd6349f6 100644 --- a/opm/autodiff/WellHelpers.hpp +++ b/opm/autodiff/WellHelpers.hpp @@ -35,7 +35,7 @@ namespace Opm { namespace wellhelpers { - + inline double rateToCompare(const std::vector& well_phase_flow_rate, const int well, @@ -147,6 +147,15 @@ namespace Opm { return dp; } + inline + double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth, + const double rho, const double gravity) { + const double dh = vfp_ref_depth - well_ref_depth; + const double dp = rho * gravity * dh; + + return dp; + } + template inline Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp new file mode 100644 index 000000000..950fca584 --- /dev/null +++ b/opm/autodiff/WellInterface.hpp @@ -0,0 +1,313 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_WELLINTERFACE_HEADER_INCLUDED +#define OPM_WELLINTERFACE_HEADER_INCLUDED + +#include + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace Opm +{ + + + template + class WellInterface + { + public: + + using WellState = WellStateFullyImplicitBlackoilDense; + + typedef BlackoilModelParameters ModelParameters; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + + static const int numEq = BlackoilIndices::numEq; + typedef double Scalar; + + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + typedef DenseAd::Evaluation Eval; + + typedef Ewoms::BlackOilPolymerModule PolymerModule; + + static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); + static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); + + /// Constructor + WellInterface(const Well* well, const int time_step, const Wells* wells); + + /// Virutal destructor + virtual ~WellInterface() {} + + /// Well name. + const std::string& name() const; + + /// Well type, INJECTOR or PRODUCER. + WellType wellType() const; + + /// Well controls + WellControls* wellControls() const; + + void setVFPProperties(const VFPProperties* vfp_properties_arg); + + virtual void init(const PhaseUsage* phase_usage_arg, + const std::vector* active_arg, + const std::vector& depth_arg, + const double gravity_arg, + const int num_cells); + + virtual void initPrimaryVariablesEvaluation() const = 0; + + /// a struct to collect information about the convergence checking + struct ConvergenceReport { + struct ProblemWell { + std::string well_name; + std::string phase_name; + }; + bool converged = true; + bool nan_residual_found = false; + std::vector nan_residual_wells; + // We consider Inf is large residual here + bool too_large_residual_found = false; + std::vector too_large_residual_wells; + + ConvergenceReport& operator+=(const ConvergenceReport& rhs) { + converged = converged && rhs.converged; + nan_residual_found = nan_residual_found || rhs.nan_residual_found; + if (rhs.nan_residual_found) { + for (const ProblemWell& well : rhs.nan_residual_wells) { + nan_residual_wells.push_back(well); + } + } + too_large_residual_found = too_large_residual_found || rhs.too_large_residual_found; + if (rhs.too_large_residual_found) { + for (const ProblemWell& well : rhs.too_large_residual_wells) { + too_large_residual_wells.push_back(well); + } + } + return *this; + } + }; + + virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg, + const ModelParameters& param) const = 0; + + virtual void solveEqAndUpdateWellState(const ModelParameters& param, + WellState& well_state) = 0; + + virtual void assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) = 0; + + void updateListEconLimited(const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const; + + void setWellEfficiencyFactor(const double efficiency_factor); + + void computeRepRadiusPerfLength(const Grid& grid, const std::map& cartesian_to_compressed); + + /// using the solution x to recover the solution xw for wells and applying + /// xw to update Well State + virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param, + WellState& well_state) const = 0; + + /// Ax = Ax - C D^-1 B x + virtual void apply(const BVector& x, BVector& Ax) const = 0; + + /// r = r - C D^-1 Rw + virtual void apply(BVector& r) const = 0; + + virtual void computeWellPotentials(const Simulator& ebosSimulator, + const WellState& well_state, + std::vector& well_potentials) const = 0; + + virtual void computeAccumWell() = 0; + + // TODO: it should come with a different name + // for MS well, the definition is different and should not use this name anymore + virtual void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw) = 0; + + virtual void updateWellStateWithTarget(const int current, + WellState& xw) const = 0; + + virtual void updateWellControl(WellState& xw, + wellhelpers::WellSwitchingLogger& logger) const = 0; + + virtual void updatePrimaryVariables(const WellState& well_state) const = 0; + + protected: + + // to indicate a invalid connection + static const int INVALIDCONNECTION = -100000; + + const Well* well_ecl_; + + const int current_step_; + + // the index of well in Wells struct + int index_of_well_; + + // well type + // INJECTOR or PRODUCER + enum WellType well_type_; + + // number of phases + int number_of_phases_; + + // component fractions for each well + // typically, it should apply to injection wells + std::vector comp_frac_; + + // controls for this well + // TODO: later will check whehter to let it stay with pointer + struct WellControls* well_controls_; + + // number of the perforations for this well + int number_of_perforations_; + + // record the index of the first perforation + // TODO: it might not be needed if we refactor WellState to be a vector + // of states of individual well. + int first_perf_; + + // well index for each perforation + std::vector well_index_; + + // TODO: it might should go to StandardWell + // depth for each perforation + std::vector perf_depth_; + + // reference depth for the BHP + double ref_depth_; + + double well_efficiency_factor_; + + // cell index for each well perforation + std::vector well_cells_; + + // saturation table nubmer for each well perforation + std::vector saturation_table_number_; + + // representative radius of the perforations, used in shear calculation + std::vector perf_rep_radius_; + + // length of the perforations, use in shear calculation + std::vector perf_length_; + + // well bore diameter + std::vector bore_diameters_; + + const PhaseUsage* phase_usage_; + + bool getAllowCrossFlow() const; + + const std::vector* active_; + + const VFPProperties* vfp_properties_; + + double gravity_; + + const std::vector& active() const; + + const PhaseUsage& phaseUsage() const; + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const; + + int flowToEbosPvIdx( const int flowPv ) const; + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + + // TODO: it is dumplicated with StandardWellsDense + int numComponents() const; + + double wsolvent() const; + + double wpolymer() const; + + bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const; + + bool wellHasTHPConstraints() const; + + // Component fractions for each phase for the well + const std::vector& compFrac() const; + + double mostStrictBhpFromBhpLimits() const; + + // a tuple type for ratio limit check. + // first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three + // values should not be used. + // second value indicates whehter there is only one connection left. + // third value indicates the indx of the worst-offending connection. + // the last value indicates the extent of the violation for the worst-offending connection, which is defined by + // the ratio of the actual value to the value of the violated limit. + using RatioCheckTuple = std::tuple; + + RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const; + + RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const; + + }; + +} + +#include "WellInterface_impl.hpp" + +#endif // OPM_WELLINTERFACE_HEADER_INCLUDED diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp new file mode 100644 index 000000000..6133291e0 --- /dev/null +++ b/opm/autodiff/WellInterface_impl.hpp @@ -0,0 +1,768 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +namespace Opm +{ + + + template + WellInterface:: + WellInterface(const Well* well, const int time_step, const Wells* wells) + : well_ecl_(well) + , current_step_(time_step) + { + if (!well) { + OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface"); + } + + if (time_step < 0) { + OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface"); + } + + if (!wells) { + OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface"); + } + + const std::string& well_name = well->name(); + + // looking for the location of the well in the wells struct + int index_well; + for (index_well = 0; index_well < wells->number_of_wells; ++index_well) { + if (well_name == std::string(wells->name[index_well])) { + break; + } + } + + // should not enter the constructor if the well does not exist in the wells struct + // here, just another assertion. + assert(index_well != wells->number_of_wells); + + index_of_well_ = index_well; + well_type_ = wells->type[index_well]; + number_of_phases_ = wells->number_of_phases; + + // copying the comp_frac + { + comp_frac_.resize(number_of_phases_); + const int index_begin = index_well * number_of_phases_; + std::copy(wells->comp_frac + index_begin, + wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() ); + } + + well_controls_ = wells->ctrls[index_well]; + + ref_depth_ = wells->depth_ref[index_well]; + + // perforations related + { + const int perf_index_begin = wells->well_connpos[index_well]; + const int perf_index_end = wells->well_connpos[index_well + 1]; + number_of_perforations_ = perf_index_end - perf_index_begin; + first_perf_ = perf_index_begin; + + well_cells_.resize(number_of_perforations_); + std::copy(wells->well_cells + perf_index_begin, + wells->well_cells + perf_index_end, + well_cells_.begin() ); + + well_index_.resize(number_of_perforations_); + std::copy(wells->WI + perf_index_begin, + wells->WI + perf_index_end, + well_index_.begin() ); + + saturation_table_number_.resize(number_of_perforations_); + std::copy(wells->sat_table_id + perf_index_begin, + wells->sat_table_id + perf_index_end, + saturation_table_number_.begin() ); + } + + well_efficiency_factor_ = 1.0; + } + + + + + + template + void + WellInterface:: + init(const PhaseUsage* phase_usage_arg, + const std::vector* active_arg, + const std::vector& /* depth_arg */, + const double gravity_arg, + const int /* num_cells */) + { + phase_usage_ = phase_usage_arg; + active_ = active_arg; + gravity_ = gravity_arg; + } + + + + + + template + void + WellInterface:: + setVFPProperties(const VFPProperties* vfp_properties_arg) + { + vfp_properties_ = vfp_properties_arg; + } + + + + + + template + const std::string& + WellInterface:: + name() const + { + return well_ecl_->name(); + } + + + + + + template + WellType + WellInterface:: + wellType() const + { + return well_type_; + } + + + + + + template + WellControls* + WellInterface:: + wellControls() const + { + return well_controls_; + } + + + + + + template + bool + WellInterface:: + getAllowCrossFlow() const + { + return well_ecl_->getAllowCrossFlow(); + } + + + + + + template + const std::vector& + WellInterface:: + active() const + { + assert(active_); + + return *active_; + } + + + + + + template + void + WellInterface:: + setWellEfficiencyFactor(const double efficiency_factor) + { + well_efficiency_factor_ = efficiency_factor; + } + + + + + + template + const PhaseUsage& + WellInterface:: + phaseUsage() const + { + assert(phase_usage_); + + return *phase_usage_; + } + + + + + + template + int + WellInterface:: + flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; + if (phaseIdx > 2 ) + return phaseIdx; + return phaseToComp[ phaseIdx ]; + } + + + + + + template + int + WellInterface:: + flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + + if (flowPv > 2 ) + return flowPv; + + return flowToEbos[ flowPv ]; + } + + + + + + template + int + WellInterface:: + flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + assert(phaseIdx < 3); + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + return flowToEbos[ phaseIdx ]; + } + + + + + + template + int + WellInterface:: + numComponents() const + { + // TODO: how about two phase polymer + if (number_of_phases_ == 2) { + return 2; + } + + int numComp = FluidSystem::numComponents; + + if (has_solvent) { + numComp++; + } + return numComp; + } + + + + + template + double + WellInterface:: + wsolvent() const + { + if (!has_solvent) { + return 0.0; + } + + WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_); + if (injection.injectorType == WellInjector::GAS) { + double solvent_fraction = well_ecl_->getSolventFraction(current_step_); + return solvent_fraction; + } + + assert(false); + return 0.0; + } + + + + + + template + double + WellInterface:: + wpolymer() const + { + if (!has_polymer) { + return 0.0; + } + + WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_); + WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_); + + if (injection.injectorType == WellInjector::WATER) { + const double polymer_injection_concentration = polymer.m_polymerConcentration; + return polymer_injection_concentration; + } + + assert(false); // TODO: find a more logical way to handle this situation + return 0.0; + } + + + + + + template + double + WellInterface:: + mostStrictBhpFromBhpLimits() const + { + double bhp; + + // initial bhp value, making the value not usable + switch( well_type_ ) { + case INJECTOR: + bhp = std::numeric_limits::max(); + break; + case PRODUCER: + bhp = -std::numeric_limits::max(); + break; + default: + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name()); + } + + // The number of the well controls/constraints + const int nwc = well_controls_get_num(well_controls_); + + for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { + // finding a BHP constraint + if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) { + // get the bhp constraint value, it should always be postive assummingly + const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index); + + switch(well_type_) { + case INJECTOR: // using the lower bhp contraint from Injectors + if (bhp_target < bhp) { + bhp = bhp_target; + } + break; + case PRODUCER: + if (bhp_target > bhp) { + bhp = bhp_target; + } + break; + default: + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name()); + } // end of switch + } + } + + return bhp; + } + + + + + template + bool + WellInterface:: + wellHasTHPConstraints() const + { + const int nwc = well_controls_get_num(well_controls_); + for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(well_controls_, ctrl_index) == THP) { + return true; + } + } + return false; + } + + + + + + template + bool + WellInterface:: + checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const + { + const Opm::PhaseUsage& pu = *phase_usage_; + const int np = number_of_phases_; + + if (econ_production_limits.onMinOilRate()) { + assert(active()[Oil]); + const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; + const double min_oil_rate = econ_production_limits.minOilRate(); + if (std::abs(oil_rate) < min_oil_rate) { + return true; + } + } + + if (econ_production_limits.onMinGasRate() ) { + assert(active()[Gas]); + const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ]; + const double min_gas_rate = econ_production_limits.minGasRate(); + if (std::abs(gas_rate) < min_gas_rate) { + return true; + } + } + + if (econ_production_limits.onMinLiquidRate() ) { + assert(active()[Oil]); + assert(active()[Water]); + const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + const double min_liquid_rate = econ_production_limits.minLiquidRate(); + if (std::abs(liquid_rate) < min_liquid_rate) { + return true; + } + } + + if (econ_production_limits.onMinReservoirFluidRate()) { + OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); + } + + return false; + } + + + + + + + template + typename WellInterface::RatioCheckTuple + WellInterface:: + checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const + { + bool water_cut_limit_violated = false; + int worst_offending_connection = INVALIDCONNECTION; + bool last_connection = false; + double violation_extent = -1.0; + + const int np = number_of_phases_; + const Opm::PhaseUsage& pu = *phase_usage_; + const int well_number = index_of_well_; + + assert(active()[Oil]); + assert(active()[Water]); + + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + double water_cut; + if (std::abs(liquid_rate) != 0.) { + water_cut = water_rate / liquid_rate; + } else { + water_cut = 0.0; + } + + const double max_water_cut_limit = econ_production_limits.maxWaterCut(); + if (water_cut > max_water_cut_limit) { + water_cut_limit_violated = true; + } + + if (water_cut_limit_violated) { + // need to handle the worst_offending_connection + const int perf_start = first_perf_; + const int perf_number = number_of_perforations_; + + std::vector water_cut_perf(perf_number); + for (int perf = 0; perf < perf_number; ++perf) { + const int i_perf = perf_start + perf; + const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; + const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; + const double liquid_perf_rate = oil_perf_rate + water_perf_rate; + if (std::abs(liquid_perf_rate) != 0.) { + water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; + } else { + water_cut_perf[perf] = 0.; + } + } + + last_connection = (perf_number == 1); + if (last_connection) { + worst_offending_connection = 0; + violation_extent = water_cut_perf[0] / max_water_cut_limit; + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + double max_water_cut_perf = 0.; + for (int perf = 0; perf < perf_number; ++perf) { + if (water_cut_perf[perf] > max_water_cut_perf) { + worst_offending_connection = perf; + max_water_cut_perf = water_cut_perf[perf]; + } + } + + assert(max_water_cut_perf != 0.); + assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); + + violation_extent = max_water_cut_perf / max_water_cut_limit; + } + + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + + + + + template + typename WellInterface::RatioCheckTuple + WellInterface:: + checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state) const + { + // TODO: not sure how to define the worst-offending connection when more than one + // ratio related limit is violated. + // The defintion used here is that we define the violation extent based on the + // ratio between the value and the corresponding limit. + // For each violated limit, we decide the worst-offending connection separately. + // Among the worst-offending connections, we use the one has the biggest violation + // extent. + + bool any_limit_violated = false; + bool last_connection = false; + int worst_offending_connection = INVALIDCONNECTION; + double violation_extent = -1.0; + + if (econ_production_limits.onMaxWaterCut()) { + const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state); + bool water_cut_violated = std::get<0>(water_cut_return); + if (water_cut_violated) { + any_limit_violated = true; + const double violation_extent_water_cut = std::get<3>(water_cut_return); + if (violation_extent_water_cut > violation_extent) { + violation_extent = violation_extent_water_cut; + worst_offending_connection = std::get<2>(water_cut_return); + last_connection = std::get<1>(water_cut_return); + } + } + } + + if (econ_production_limits.onMaxGasOilRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxWaterGasRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxGasLiquidRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); + } + + if (any_limit_violated) { + assert(worst_offending_connection >=0); + assert(violation_extent > 1.); + } + + return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + + + + + template + void + WellInterface:: + updateListEconLimited(const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const + { + // economic limits only apply for production wells. + if (wellType() != PRODUCER) { + return; + } + + // flag to check if the mim oil/gas rate limit is violated + bool rate_limit_violated = false; + const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_); + + // if no limit is effective here, then continue to the next well + if ( !econ_production_limits.onAnyEffectiveLimit() ) { + return; + } + + const std::string well_name = name(); + + // for the moment, we only handle rate limits, not handling potential limits + // the potential limits should not be difficult to add + const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); + if (quantity_limit == WellEcon::POTN) { + const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") + + std::string("All the limits will be evaluated based on RATE. "); + OpmLog::warning("NOT_SUPPORTING_POTN", msg); + } + + if (econ_production_limits.onAnyRateLimit()) { + rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state); + } + + if (rate_limit_violated) { + if (econ_production_limits.endRun()) { + const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") + + std::string("the program will keep running after ") + well_name + std::string(" is closed"); + OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); + } + + if (econ_production_limits.validFollowonWell()) { + OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); + } + + if (well_ecl_->getAutomaticShutIn()) { + list_econ_limited.addShutWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); + OpmLog::info(msg); + } else { + list_econ_limited.addStoppedWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); + OpmLog::info(msg); + } + // the well is closed, not need to check other limits + return; + } + + // checking for ratio related limits, mostly all kinds of ratio. + bool ratio_limits_violated = false; + RatioCheckTuple ratio_check_return; + + if (econ_production_limits.onAnyRatioLimit()) { + ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state); + ratio_limits_violated = std::get<0>(ratio_check_return); + } + + if (ratio_limits_violated) { + const bool last_connection = std::get<1>(ratio_check_return); + const int worst_offending_connection = std::get<2>(ratio_check_return); + + assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_)); + + const int cell_worst_offending_connection = well_cells_[worst_offending_connection]; + list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); + const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") + + well_name + std::string(" will be closed due to economic limit"); + OpmLog::info(msg); + + if (last_connection) { + list_econ_limited.addShutWell(well_name); + const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); + OpmLog::info(msg2); + } + } + } + + + + + + template + void + WellInterface:: + computeRepRadiusPerfLength(const Grid& grid, + const std::map& cartesian_to_compressed) + { + const int* cart_dims = Opm::UgGridHelpers::cartDims(grid); + auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid); + auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid); + + const int nperf = number_of_perforations_; + + perf_rep_radius_.clear(); + perf_length_.clear(); + bore_diameters_.clear(); + + perf_rep_radius_.resize(nperf); + perf_length_.resize(nperf); + bore_diameters_.resize(nperf); + + // COMPDAT handling + const auto& completionSet = well_ecl_->getCompletions(current_step_); + for (size_t c=0; c::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << name() << ')'); + } + const int cell = cgit->second; + + { + double radius = 0.5*completion.getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + + const std::array cubical = + WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell); + + WellCompletion::DirectionEnum direction = completion.getDirection(); + + double re; // area equivalent radius of the grid block + double perf_length; // the length of the well perforation + + switch (direction) { + case Opm::WellCompletion::DirectionEnum::X: + re = std::sqrt(cubical[1] * cubical[2] / M_PI); + perf_length = cubical[0]; + break; + case Opm::WellCompletion::DirectionEnum::Y: + re = std::sqrt(cubical[0] * cubical[2] / M_PI); + perf_length = cubical[1]; + break; + case Opm::WellCompletion::DirectionEnum::Z: + re = std::sqrt(cubical[0] * cubical[1] / M_PI); + perf_length = cubical[2]; + break; + default: + OPM_THROW(std::runtime_error, " Dirtecion of well is not supported "); + } + + const double repR = std::sqrt(re * radius); + perf_rep_radius_.push_back(repR); + perf_length_.push_back(perf_length); + bore_diameters_.push_back(2. * radius); + } + } + } + } + +} diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index 48de2d124..d7b6d499f 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -108,101 +108,6 @@ namespace Opm } } } - - - // TODO: the reason to keep this is to avoid getting defaulted value BHP - // some facilities needed from opm-parser or opm-core - // It is a little tricky, since sometimes before applying group control, the only - // available constraints in the well_controls is the defaulted BHP value, and it - // is really not desirable to use this value to enter the Newton iterations. - setWellSolutions(pu); - } - - - - /// Set wellSolutions() based on the base class members. - void setWellSolutions(const PhaseUsage& pu) - { - // Set nw and np, or return if no wells. - if (wells_.get() == nullptr) { - return; - } - const int nw = wells_->number_of_wells; - if (nw == 0) { - return; - } - - - const int np = wells_->number_of_phases; - const int numComp = pu.has_solvent? np+1:np; - well_solutions_.clear(); - well_solutions_.resize(nw * numComp, 0.0); - std::vector g = {1.0,1.0,0.01}; - for (int w = 0; w < nw; ++w) { - WellControls* wc = wells_->ctrls[w]; - - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - const int current = currentControls()[w]; - well_controls_set_current( wc, current); - const WellType& well_type = wells_->type[w]; - - switch (well_controls_iget_type(wc, current)) { - case THP: // Intentional fall-through - case BHP: - if (well_type == INJECTOR) { - for (int p = 0; p < np; ++p) { - well_solutions_[w] += wellRates()[np*w + p] * wells_->comp_frac[np*w + p]; - } - } else { - for (int p = 0; p < np; ++p) { - well_solutions_[w] += g[p] * wellRates()[np*w + p]; - } - } - break; - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - wellSolutions()[w] = bhp()[w]; - break; - } - - double total_rates = 0.0; - for (int p = 0; p < np; ++p) { - total_rates += g[p] * wellRates()[np*w + p]; - } - - const int waterpos = pu.phase_pos[Water]; - const int gaspos = pu.phase_pos[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] ) { - wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates; - } - if( pu.phase_used[Gas] ) { - wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ; - } - if( pu.has_solvent) { - wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates; - } - - - - } else { - if( pu.phase_used[Water] ) { - wellSolutions()[nw + w] = wells_->comp_frac[np*w + waterpos]; - } - if( pu.phase_used[Gas] ) { - wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos]; - } - if (pu.has_solvent) { - wellSolutions()[3*nw + w] = 0; - } - - } - } } @@ -212,11 +117,6 @@ namespace Opm init(wells, state.pressure(), dummy_state, pu) ; } - - /// One rate per phase and well connection. - std::vector& wellSolutions() { return well_solutions_; } - const std::vector& wellSolutions() const { return well_solutions_; } - /// One rate pr well connection. std::vector& perfRateSolvent() { return perfRateSolvent_; } const std::vector& perfRateSolvent() const { return perfRateSolvent_; } @@ -252,7 +152,6 @@ namespace Opm private: - std::vector well_solutions_; std::vector perfRateSolvent_; }; diff --git a/tests/TESTWELLMODEL.DATA b/tests/TESTWELLMODEL.DATA new file mode 100644 index 000000000..796b4ff68 --- /dev/null +++ b/tests/TESTWELLMODEL.DATA @@ -0,0 +1,58 @@ +RUNSPEC + +DIMENS + 5 5 4 / +OIL +GAS +WATER + +GRID + +DX + 100*100. / + +DY + 100*50. / + +DZ + 100*10. / + +TOPS + 25*2500 / + +PORO + 100*0.3 / + +PERMX + 100*10. / + +PERMY + 100*20. / + +PERMZ + 100*1. / + +SCHEDULE + +WELSPECS + 'PROD1' 'P' 5 5 2525 'OIL' / + 'INJE1' 'I' 1 1 2505 'GAS' / +/ + +COMPDAT + 'PROD1' 5 5 3 4 'OPEN' 2* 0.15 / + 'INJE1' 1 1 1 4 'OPEN' 2* 0.15 / +/ + +WCONINJE + 'INJE1' 'WATER' 'OPEN' 'RATE' 1000. / +/ + +WCONPROD + 'PROD1' 'OPEN' 'GRAT' 2* 50000. / +/ + +TSTEP +10 / + +END diff --git a/tests/test_multisegmentwells.cpp b/tests/test_multisegmentwells.cpp index 82a00c39b..403aad311 100644 --- a/tests/test_multisegmentwells.cpp +++ b/tests/test_multisegmentwells.cpp @@ -83,15 +83,6 @@ struct SetupMSW { std::unique_ptr grid_init(new GridInit(ecl_state, porv)); const Grid& grid = grid_init->grid(); - // Create material law manager. - std::vector compressed_to_cartesianIdx; - Opm::createGlobalCellArray(grid, compressed_to_cartesianIdx); - - std::shared_ptr material_law_manager(new MaterialLawManager()); - material_law_manager->initFromDeck(deck, ecl_state, compressed_to_cartesianIdx); - - std::unique_ptr fluidprops(new FluidProps(deck, ecl_state, material_law_manager, grid)); - const size_t current_timestep = 0; // dummy_dynamic_list_econ_lmited @@ -115,7 +106,7 @@ struct SetupMSW { std::unordered_set()); const Wells* wells = wells_manager.c_wells(); - const auto wells_ecl = ecl_state.getSchedule().getWells(current_timestep); + const auto& wells_ecl = ecl_state.getSchedule().getWells(current_timestep); ms_wells.reset(new Opm::MultisegmentWells(wells, &(wells_manager.wellCollection()), wells_ecl, current_timestep)); }; diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp new file mode 100644 index 000000000..c954e8e0a --- /dev/null +++ b/tests/test_wellmodel.cpp @@ -0,0 +1,194 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE WellModelTest + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +// maybe should just include BlackoilModelEbos.hpp +namespace Ewoms { + namespace Properties { + NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); + } +} + +using StandardWell = Opm::StandardWell; + +struct SetupTest { + + using Grid = UnstructuredGrid; + using GridInit = Opm::GridInit; + + SetupTest () + { + Opm::ParseContext parse_context; + Opm::Parser parser; + auto deck = parser.parseFile("TESTWELLMODEL.DATA", parse_context); + ecl_state = std::make_unique(deck , parse_context); + + // Create grid. + const std::vector& porv = + ecl_state->get3DProperties().getDoubleGridProperty("PORV").getData(); + + std::unique_ptr grid_init(new GridInit(*ecl_state, porv)); + const Grid& grid = grid_init->grid(); + + // Create material law manager. + std::vector compressed_to_cartesianIdx; + Opm::createGlobalCellArray(grid, compressed_to_cartesianIdx); + + // dummy_dynamic_list_econ_lmited + const Opm::DynamicListEconLimited dummy_dynamic_list; + + current_timestep = 0; + + // Create wells. + wells_manager = std::make_unique (*ecl_state, + current_timestep, + Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::cartDims(grid), + Opm::UgGridHelpers::dimensions(grid), + Opm::UgGridHelpers::cell2Faces(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + dummy_dynamic_list, + false, + std::unordered_set()); + + }; + + std::unique_ptr wells_manager; + std::unique_ptr ecl_state; + int current_timestep; +}; + + +BOOST_AUTO_TEST_CASE(TestStandardWellInput) { + SetupTest setup_test; + const Wells* wells = setup_test.wells_manager->c_wells(); + const auto& wells_ecl = setup_test.ecl_state->getSchedule().getWells(setup_test.current_timestep); + BOOST_CHECK_EQUAL( wells_ecl.size(), 2); + const Opm::Well* well = wells_ecl[1]; + BOOST_CHECK_THROW( StandardWell( well, -1, wells), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( well, 4, nullptr), std::invalid_argument); +} + + +BOOST_AUTO_TEST_CASE(TestBehavoir) { + SetupTest setup_test; + const Wells* wells_struct = setup_test.wells_manager->c_wells(); + const auto& wells_ecl = setup_test.ecl_state->getSchedule().getWells(setup_test.current_timestep); + const int current_timestep = setup_test.current_timestep; + std::vector > wells; + + { + const int nw = wells_struct ? (wells_struct->number_of_wells) : 0; + + for (int w = 0; w < nw; ++w) { + const std::string well_name(wells_struct->name[w]); + + size_t index_well = 0; + for (; index_well < wells_ecl.size(); ++index_well) { + if (well_name == wells_ecl[index_well]->name()) { + break; + } + } + // we should always be able to find the well in wells_ecl + BOOST_CHECK(index_well != wells_ecl.size()); + + wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct) ); + } + } + + // first well, it is a production well from the deck + { + const auto& well = wells[0]; + BOOST_CHECK_EQUAL(well->name(), "PROD1"); + BOOST_CHECK(well->wellType() == PRODUCER); + BOOST_CHECK(well->numEq == 3); + BOOST_CHECK(well->numWellEq == 3); + const auto& wc = well->wellControls(); + const int ctrl_num = well_controls_get_num(wc); + BOOST_CHECK(ctrl_num > 0); + const auto& control = well_controls_get_current(wc); + BOOST_CHECK(control >= 0); + // GAS RATE CONTROL + const auto& distr = well_controls_iget_distr(wc, control); + BOOST_CHECK(distr[0] == 0.); + BOOST_CHECK(distr[1] == 0.); + BOOST_CHECK(distr[2] == 1.); + } + + // second well, it is the injection well from the deck + { + const auto& well = wells[1]; + BOOST_CHECK_EQUAL(well->name(), "INJE1"); + BOOST_CHECK(well->wellType() == INJECTOR); + BOOST_CHECK(well->numEq == 3); + BOOST_CHECK(well->numWellEq == 3); + const auto& wc = well->wellControls(); + const int ctrl_num = well_controls_get_num(wc); + BOOST_CHECK(ctrl_num > 0); + const auto& control = well_controls_get_current(wc); + BOOST_CHECK(control >= 0); + // WATER RATE CONTROL + const auto& distr = well_controls_iget_distr(wc, control); + BOOST_CHECK(distr[0] == 1.); + BOOST_CHECK(distr[1] == 0.); + BOOST_CHECK(distr[2] == 0.); + } +}