diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index e162b1de5..e73e6fe8e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -250,7 +250,11 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellInterface_impl.hpp opm/autodiff/StandardWell.hpp opm/autodiff/StandardWell_impl.hpp - opm/autodiff/StandardWellsDense.hpp + opm/autodiff/MultisegmentWell.hpp + opm/autodiff/MultisegmentWell_impl.hpp + opm/autodiff/MSWellHelpers.hpp + opm/autodiff/BlackoilWellModel.hpp + opm/autodiff/BlackoilWellModel_impl.hpp opm/autodiff/StandardWellsSolvent.hpp opm/autodiff/StandardWellsSolvent_impl.hpp opm/autodiff/MissingFeatures.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index cc2b43802..237b624b8 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -28,7 +28,7 @@ #include #include -#include +#include #include #include #include @@ -154,7 +154,7 @@ namespace Opm { /// \param[in] terminal_output request output to cout/cerr BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, - StandardWellsDense& well_model, + BlackoilWellModel& well_model, RateConverterType& rate_converter, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output @@ -391,7 +391,7 @@ namespace Opm { } catch ( const Dune::FMatrixError& e ) { - OPM_THROW(Opm::NumericalProblem,"Well equation did not converge"); + OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations"); } return report; @@ -425,12 +425,12 @@ namespace Opm { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx]; } - + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx]; - oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx]; + oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx]; } - + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew; } @@ -446,16 +446,16 @@ namespace Opm { saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; } - + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; } - + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; } - + Scalar tmp = pressureNew - pressureOld; resultDelta += tmp*tmp; resultDenom += pressureNew*pressureNew; @@ -503,14 +503,14 @@ namespace Opm { // Solve system. if( isParallel() ) { - typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense, true > Operator; + typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel, true > Operator; Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() ); assert( opA.comm() ); istlSolver().solve( opA, x, ebosResid, *(opA.comm()) ); } else { - typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense, false > Operator; + typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel, false > Operator; Operator opA(ebosJac, well_model_); istlSolver().solve( opA, x, ebosResid ); } @@ -680,7 +680,7 @@ namespace Opm { maxVal = std::max(std::abs(dsw),maxVal); maxVal = std::max(std::abs(dsg),maxVal); maxVal = std::max(std::abs(dso),maxVal); - maxVal = std::max(std::abs(dss),maxVal); + maxVal = std::max(std::abs(dss),maxVal); double satScaleFactor = 1.0; if (maxVal > dsMax()) { @@ -905,7 +905,6 @@ namespace Opm { bool converged_MB = true; bool converged_CNV = true; - bool converged_Well = true; // Finish computation for ( int compIdx = 0; compIdx < numComp; ++compIdx ) { @@ -913,13 +912,12 @@ namespace Opm { mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum; converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb); converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv); - // Well flux convergence is only for fluid phases, not other materials - // in our current implementation. - converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg); residual_norms.push_back(CNV[compIdx]); } + const bool converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg); + bool converged = converged_MB && converged_Well; // do not care about the cell based residual in the last two Newton @@ -1502,7 +1500,7 @@ namespace Opm { SimulatorReport failureReport_; // Well Model - StandardWellsDense& well_model_; + BlackoilWellModel& well_model_; /// \brief Whether we print something to std::cout bool terminal_output_; @@ -1519,9 +1517,10 @@ namespace Opm { public: /// return the StandardWells object - StandardWellsDense& + BlackoilWellModel& wellModel() { return well_model_; } - const StandardWellsDense& + + const BlackoilWellModel& wellModel() const { return well_model_; } int numWells() const { return well_model_.numWells(); } diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 9808556a7..a910e9c15 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -50,12 +50,20 @@ namespace Opm tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_); tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ ); tolerance_well_control_ = param.getDefault("tolerance_well_control", tolerance_well_control_); + max_welleq_iter_ = param.getDefault("max_welleq_iter", max_welleq_iter_); + if (use_multisegment_well_) { + tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_); + max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_); + use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_); + max_inner_iter_ms_wells_ = param.getDefault("max_inner_iter_ms_wells", max_inner_iter_ms_wells_); + } maxSinglePrecisionTimeStep_ = unit::convert::from( param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_); update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_); use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_); + use_multisegment_well_ = param.getDefault("use_multisegment_well", use_multisegment_well_); deck_file_name_ = param.template get("deck_filename"); } @@ -75,10 +83,16 @@ namespace Opm tolerance_cnv_ = 1.0e-2; tolerance_wells_ = 1.0e-4; tolerance_well_control_ = 1.0e-7; + tolerance_pressure_ms_wells_ = unit::convert::from(0.01, unit::barsa); // 0.01 bar + max_welleq_iter_ = 15; + max_pressure_change_ms_wells_ = unit::convert::from(2.0, unit::barsa); // 2.0 bar + use_inner_iterations_ms_wells_ = true; + max_inner_iter_ms_wells_ = 10; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; update_equations_scaling_ = false; use_update_stabilization_ = true; + use_multisegment_well_ = false; } diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index cfb36728b..b14ab1135 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -51,6 +51,19 @@ namespace Opm /// Tolerance for the well control equations // TODO: it might need to distinguish between rate control and pressure control later double tolerance_well_control_; + /// Tolerance for the pressure equations for multisegment wells + double tolerance_pressure_ms_wells_; + /// Maximum pressure change over an iteratio for ms wells + double max_pressure_change_ms_wells_; + + /// Whether to use inner iterations for ms wells + bool use_inner_iterations_ms_wells_; + + /// Maximum inner iteration number for ms wells + int max_inner_iter_ms_wells_; + + /// Maximum iteration number of the well equation solution + int max_welleq_iter_; /// Tolerance for time step in seconds where single precision can be used /// for solving for the Jacobian @@ -68,7 +81,14 @@ namespace Opm /// Try to detect oscillation or stagnation. bool use_update_stabilization_; - // The file name of the deck + /// Whether to use MultisegmentWell to handle multisegment wells + /// it is something temporary before the multisegment well model is considered to be + /// well developed and tested. + /// if it is false, we will handle multisegment wells as standard wells, which will be + /// the default behavoir for the moment. Later, we might set it to be true by default if necessary + bool use_multisegment_well_; + + /// The file name of the deck std::string deck_file_name_; /// Construct from user parameters or defaults. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/BlackoilWellModel.hpp similarity index 82% rename from opm/autodiff/StandardWellsDense.hpp rename to opm/autodiff/BlackoilWellModel.hpp index 078706afb..e152cbf36 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -21,8 +21,8 @@ */ -#ifndef OPM_STANDARDWELLSDENSE_HEADER_INCLUDED -#define OPM_STANDARDWELLSDENSE_HEADER_INCLUDED +#ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED #include @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -60,9 +61,9 @@ namespace Opm { - /// Class for handling the standard well model. + /// Class for handling the blackoil well model. template - class StandardWellsDense { + class BlackoilWellModel { public: // --------- Types --------- typedef WellStateFullyImplicitBlackoil WellState; @@ -90,14 +91,14 @@ namespace Opm { SurfaceToReservoirVoidage >; // --------- Public methods --------- - StandardWellsDense(const Wells* wells_arg, - WellCollection* well_collection, - const std::vector< const Well* >& wells_ecl, - const ModelParameters& param, - const RateConverterType& rate_converter, - const bool terminal_output, - const int current_index, - std::vector& pvt_region_idx); + BlackoilWellModel(const Wells* wells_arg, + WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, + const ModelParameters& param, + const RateConverterType& rate_converter, + const bool terminal_output, + const int current_index, + const std::vector& pvt_region_idx); void init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, @@ -137,7 +138,7 @@ namespace Opm { /// return true if wells are available on this process bool localWellsActive() const; - bool getWellConvergence(Simulator& ebosSimulator, + bool getWellConvergence(const Simulator& ebosSimulator, const std::vector& B_avg) const; /// upate the dynamic lists related to economic limits @@ -162,6 +163,8 @@ namespace Opm { const int number_of_phases_; + const ModelParameters& param_; + using WellInterfacePtr = std::unique_ptr >; // a vector of all the wells. // eventually, the wells_ above should be gone. @@ -174,12 +177,13 @@ namespace Opm { // create the well container static std::vector createWellContainer(const Wells* wells, const std::vector& wells_ecl, - const int time_step); + const bool use_multisegment_well, + const int time_step, + const ModelParameters& param); // Well collection is used to enforce the group control WellCollection* well_collection_; - ModelParameters param_; bool terminal_output_; bool has_solvent_; bool has_polymer_; @@ -188,7 +192,7 @@ namespace Opm { PhaseUsage phase_usage_; std::vector active_; const RateConverterType& rate_converter_; - std::vector pvt_region_idx_; + const std::vector& pvt_region_idx_; // the number of the cells in the local grid int number_of_cells_; @@ -210,7 +214,7 @@ namespace Opm { void computeRepRadiusPerfLength(const Grid& grid); - void computeAverageFormationFactor(Simulator& ebosSimulator, + void computeAverageFormationFactor(const Simulator& ebosSimulator, std::vector& B_avg) const; void applyVREPGroupControl(WellState& well_state) const; @@ -229,15 +233,18 @@ namespace Opm { void calculateEfficiencyFactors(); - void computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw) const; + // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation() + // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions + // twice at the beginning of the time step + /// Calculating the explict quantities used in the well calculation. By explicit, we mean they are cacluated + /// at the beginning of the time step and no derivatives are included in these quantities + void calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& xw) const; SimulatorReport solveWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state) const; - void computeAccumWells() const; - void initPrimaryVariablesEvaluation() const; // The number of components in the model. @@ -268,11 +275,15 @@ namespace Opm { // 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); + WellState& well_state) const; + + void prepareGroupControl(const Simulator& ebos_simulator, + WellState& well_state) const; + }; } // namespace Opm -#include "StandardWellsDense_impl.hpp" +#include "BlackoilWellModel_impl.hpp" #endif diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp similarity index 87% rename from opm/autodiff/StandardWellsDense_impl.hpp rename to opm/autodiff/BlackoilWellModel_impl.hpp index 19e324a67..b3721055e 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -4,23 +4,23 @@ namespace Opm { template - StandardWellsDense:: - StandardWellsDense(const Wells* wells_arg, - WellCollection* well_collection, - const std::vector< const Well* >& wells_ecl, - const ModelParameters& param, - const RateConverterType& rate_converter, - const bool terminal_output, - const int current_timeIdx, - std::vector& pvt_region_idx) + BlackoilWellModel:: + BlackoilWellModel(const Wells* wells_arg, + WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, + const ModelParameters& param, + const RateConverterType& rate_converter, + const bool terminal_output, + const int current_timeIdx, + const std::vector& pvt_region_idx) : 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) + , well_container_(createWellContainer(wells_arg, wells_ecl, param.use_multisegment_well_, current_timeIdx, param) ) + , well_collection_(well_collection) , terminal_output_(terminal_output) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) @@ -36,7 +36,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, const double gravity_arg, @@ -87,7 +87,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: setVFPProperties(const VFPProperties* vfp_properties_arg) { for (auto& well : well_container_) { @@ -102,7 +102,7 @@ namespace Opm { template int - StandardWellsDense:: + BlackoilWellModel:: numWells() const { return number_of_wells_; @@ -113,11 +113,13 @@ namespace Opm { template - std::vector::WellInterfacePtr > - StandardWellsDense:: + std::vector::WellInterfacePtr > + BlackoilWellModel:: createWellContainer(const Wells* wells, const std::vector< const Well* >& wells_ecl, - const int time_step) + const bool use_multisegment_well, + const int time_step, + const ModelParameters& param) { std::vector well_container; @@ -146,13 +148,12 @@ namespace Opm { } const Well* well_ecl = wells_ecl[index_well]; - // TODO: stopping throwing when encoutnering MS wells for now. - /* 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) ); + if ( !well_ecl->isMultiSegment(time_step) || !use_multisegment_well) { + well_container.emplace_back(new StandardWell(well_ecl, time_step, wells, param) ); + } else { + well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells, param) ); + } } } return well_container; @@ -164,29 +165,27 @@ namespace Opm { template SimulatorReport - StandardWellsDense:: + BlackoilWellModel:: assemble(Simulator& ebosSimulator, const int iterationIdx, const double dt, WellState& well_state) { - - if (iterationIdx == 0) { - prepareTimeStep(ebosSimulator, well_state); - } - SimulatorReport report; if ( ! wellsActive() ) { return report; } + if (iterationIdx == 0) { + prepareTimeStep(ebosSimulator, well_state); + } + updateWellControls(well_state); // Set the well primary variables based on the value of well solutions initPrimaryVariablesEvaluation(); if (iterationIdx == 0) { - computeWellConnectionPressures(ebosSimulator, well_state); - computeAccumWells(); + calculateExplicitQuantities(ebosSimulator, well_state); } if (param_.solve_welleq_initially_ && iterationIdx == 0) { @@ -205,7 +204,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: assembleWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state, @@ -224,7 +223,7 @@ namespace Opm { // r = r - duneC_^T * invDuneD_ * resWell_ template void - StandardWellsDense:: + BlackoilWellModel:: apply( BVector& r) const { if ( ! localWellsActive() ) { @@ -243,7 +242,7 @@ namespace Opm { // Ax = A x - C D^-1 B x template void - StandardWellsDense:: + BlackoilWellModel:: apply(const BVector& x, BVector& Ax) const { // TODO: do we still need localWellsActive()? @@ -263,7 +262,7 @@ namespace Opm { // Ax = Ax - alpha * C D^-1 B x template void - StandardWellsDense:: + BlackoilWellModel:: applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const { if ( ! localWellsActive() ) { @@ -287,11 +286,11 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const { for (auto& well : well_container_) { - well->recoverWellSolutionAndUpdateWellState(x, param_, well_state); + well->recoverWellSolutionAndUpdateWellState(x, well_state); } } @@ -301,7 +300,7 @@ namespace Opm { template int - StandardWellsDense:: + BlackoilWellModel:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { const auto& pu = phase_usage_; @@ -323,7 +322,7 @@ namespace Opm { template int - StandardWellsDense:: + BlackoilWellModel:: numPhases() const { return number_of_phases_; @@ -335,7 +334,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: resetWellControlFromState(const WellState& xw) const { const int nw = wells_->number_of_wells; @@ -351,7 +350,7 @@ namespace Opm { template bool - StandardWellsDense:: + BlackoilWellModel:: wellsActive() const { return wells_active_; @@ -363,7 +362,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: setWellsActive(const bool wells_active) { wells_active_ = wells_active; @@ -375,7 +374,7 @@ namespace Opm { template bool - StandardWellsDense:: + BlackoilWellModel:: localWellsActive() const { return number_of_wells_ > 0; @@ -387,7 +386,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: initPrimaryVariablesEvaluation() const { for (auto& well : well_container_) { @@ -399,23 +398,9 @@ namespace Opm { - template - void - StandardWellsDense:: - computeAccumWells() const - { - for (auto& well : well_container_) { - well->computeAccumWell(); - } - } - - - - - template SimulatorReport - StandardWellsDense:: + BlackoilWellModel:: solveWellEq(Simulator& ebosSimulator, const double dt, WellState& well_state) const @@ -427,6 +412,8 @@ namespace Opm { std::vector< Scalar > B_avg( numComp, Scalar() ); computeAverageFormationFactor(ebosSimulator, B_avg); + const int max_iter = param_.max_welleq_iter_; + int it = 0; bool converged; do { @@ -447,7 +434,7 @@ namespace Opm { if( localWellsActive() ) { for (auto& well : well_container_) { - well->solveEqAndUpdateWellState(param_, well_state); + well->solveEqAndUpdateWellState(well_state); } } // updateWellControls uses communication @@ -458,7 +445,7 @@ namespace Opm { updateWellControls(well_state); initPrimaryVariablesEvaluation(); } - } while (it < 15); + } while (it < max_iter); if (converged) { if ( terminal_output_ ) { @@ -490,14 +477,14 @@ namespace Opm { template bool - StandardWellsDense:: - getWellConvergence(Simulator& ebosSimulator, + BlackoilWellModel:: + getWellConvergence(const Simulator& ebosSimulator, const std::vector& B_avg) const { ConvergenceReport report; for (const auto& well : well_container_) { - report += well->getWellConvergence(ebosSimulator, B_avg, param_); + report += well->getWellConvergence(B_avg); } // checking NaN residuals @@ -549,14 +536,12 @@ namespace Opm { template void - StandardWellsDense:: - computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw) const + BlackoilWellModel:: + calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& xw) const { - if( ! localWellsActive() ) return ; - for (auto& well : well_container_) { - well->computeWellConnectionPressures(ebosSimulator, xw); + well->calculateExplicitQuantities(ebosSimulator, xw); } } @@ -566,7 +551,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: updateWellControls(WellState& xw) const { // Even if there no wells active locally, we cannot @@ -590,7 +575,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: updateListEconLimited(const Schedule& schedule, const int current_step, const Wells* wells_struct, @@ -608,18 +593,11 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: computeWellPotentials(const Simulator& ebosSimulator, 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 = number_of_wells_; const int np = number_of_phases_; @@ -642,18 +620,50 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: prepareTimeStep(const Simulator& ebos_simulator, - WellState& well_state) + WellState& well_state) const { - 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 - // we need to synchronize these two. - resetWellControlFromState(well_state); + // after restarting, the well_controls can be modified while + // the well_state still uses the old control index + // we need to synchronize these two. + // keep in mind that we set the control index of well_state to be the same with + // with the wellControls from the deck when we create well_state at the beginning of the report step + resetWellControlFromState(well_state); - if (wellCollection()->groupControlActive()) { + // process group control related + prepareGroupControl(ebos_simulator, well_state); + + // since the controls are all updated, we should update well_state accordingly + for (int w = 0; w < number_of_wells_; ++w) { + WellControls* wc = well_container_[w]->wellControls(); + const int control = well_controls_get_current(wc); + well_state.currentControls()[w] = control; + // 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 + if (well_state.isNewWell(w) ) { + well_state.setNewWell(w, false); + } + } // end of for (int w = 0; w < nw; ++w) + } + + + + + + template + void + BlackoilWellModel:: + prepareGroupControl(const Simulator& ebos_simulator, + WellState& well_state) const + { + // group control related processing + if (well_collection_->groupControlActive()) { + for (int w = 0; w < number_of_wells_; ++w) { WellControls* wc = well_container_[w]->wellControls(); WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name()); @@ -680,9 +690,7 @@ namespace Opm { well_node.setIndividualControl(true); } } - } - if (well_collection_->groupControlActive()) { if (well_collection_->requireWellPotentials()) { // calculate the well potentials @@ -703,22 +711,6 @@ namespace Opm { wellCollection()->updateWellTargets(well_state.wellRates()); } } - - // since the controls are all updated, we should update well_state accordingly - for (int w = 0; w < nw; ++w) { - WellControls* wc = well_container_[w]->wellControls(); - const int control = well_controls_get_current(wc); - well_state.currentControls()[w] = control; - // 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 - if (well_state.isNewWell(w) ) { - well_state.setNewWell(w, false); - } - } // end of for (int w = 0; w < nw; ++w) } @@ -727,7 +719,7 @@ namespace Opm { template WellCollection* - StandardWellsDense:: + BlackoilWellModel:: wellCollection() const { return well_collection_; @@ -739,7 +731,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: calculateEfficiencyFactors() { if ( !localWellsActive() ) { @@ -764,7 +756,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: computeWellVoidageRates(const WellState& well_state, std::vector& well_voidage_rates, std::vector& voidage_conversion_coeffs) const @@ -826,7 +818,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: applyVREPGroupControl(WellState& well_state) const { if ( wellCollection()->havingVREPGroups() ) { @@ -854,7 +846,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: updateGroupControls(WellState& well_state) const { @@ -897,7 +889,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const { if (global_cell) { @@ -915,7 +907,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: computeRepRadiusPerfLength(const Grid& grid) { // TODO, the function does not work for parallel running @@ -937,8 +929,8 @@ namespace Opm { template void - StandardWellsDense:: - computeAverageFormationFactor(Simulator& ebosSimulator, + BlackoilWellModel:: + computeAverageFormationFactor(const Simulator& ebosSimulator, std::vector& B_avg) const { const int np = numPhases(); @@ -984,7 +976,7 @@ namespace Opm { template void - StandardWellsDense:: + BlackoilWellModel:: updatePrimaryVariables(const WellState& well_state) const { for (const auto& well : well_container_) { diff --git a/opm/autodiff/Compat.hpp b/opm/autodiff/Compat.hpp index 17c4e79c3..2781a7a95 100644 --- a/opm/autodiff/Compat.hpp +++ b/opm/autodiff/Compat.hpp @@ -31,7 +31,6 @@ namespace Opm { // Forward declarations class SimulationDataContainer; class WellStateFullyImplicitBlackoil; - class WellStateFullyImplicitBlackoilDense; /// Extract single data vector from striped data. /// \return u such that u[i] = v[offset + i*stride]. @@ -71,12 +70,6 @@ namespace Opm { PhaseUsage phases, WellStateFullyImplicitBlackoil& state ); - /// As the WellStateFullyImplicitBlackoil overload, but also sets - /// the wellSolution field from the values of the other fields. - void wellsToState( const data::Wells& wells, - PhaseUsage phases, - WellStateFullyImplicitBlackoilDense& state ); - } #endif //OPM_SIMULATORS_COMPAT_HPP diff --git a/opm/autodiff/MSWellHelpers.hpp b/opm/autodiff/MSWellHelpers.hpp new file mode 100644 index 000000000..4af83165a --- /dev/null +++ b/opm/autodiff/MSWellHelpers.hpp @@ -0,0 +1,195 @@ +/* + 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_MSWELLHELPERS_HEADER_INCLUDED +#define OPM_MSWELLHELPERS_HEADER_INCLUDED + +#include +#include +#include +#include + +namespace Opm { + +namespace mswellhelpers +{ + + // obtain y = D^-1 * x with a direct solver + template + VectorType + invDXDirect(const MatrixType& D, VectorType x) + { + VectorType y(x.size()); + y = 0.; + + Dune::UMFPack linsolver(D, 0); + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult res; + + // Solve + linsolver.apply(y, x, res); + + // Checking if there is any inf or nan in y + // it will be the solution before we find a way to catch the singularity of the matrix + for (size_t i_block = 0; i_block < y.size(); ++i_block) { + for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) { + if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) { + OPM_THROW(Opm::NumericalProblem, "nan or inf value found in invDXDirect due to singular matrix"); + } + } + } + + return y; + } + + + + + + // obtain y = D^-1 * x with a BICSSTAB iterative solver + template + VectorType + invDX(const MatrixType& D, VectorType x) + { + // the function will change the value of x, so we should not use reference of x here. + + // TODO: store some of the following information to avoid to call it again and again for + // efficiency improvement. + // Bassically, only the solve / apply step is different. + + VectorType y(x.size()); + y = 0.; + + Dune::MatrixAdapter linearOperator(D); + + // Sequential incomplete LU decomposition as the preconditioner + Dune::SeqILU0 preconditioner(D, 1.0); + // Dune::SeqILUn preconditioner(D, 1, 0.92); + // Dune::SeqGS preconditioner(D, 1, 1); + // Dune::SeqJac preconditioner(D, 1, 1); + + // Preconditioned BICGSTAB solver + Dune::BiCGSTABSolver linsolver(linearOperator, + preconditioner, + 1.e-8, // desired residual reduction factor + 250, // maximum number of iterations + 0); // verbosity of the solver */ + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult res; + + // Solve + linsolver.apply(y, x, res); + + if ( !res.converged ) { + OPM_THROW(Opm::NumericalProblem, "the invDX does not get converged! "); + } + + return y; + } + + + + + + inline double haalandFormular(const double re, const double diameter, const double roughness) + { + const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) ); + + // sqrt(1/f) should be non-positive + assert(value >= 0.0); + + return 1. / (value * value); + } + + + + + + inline double calculateFrictionFactor(const double area, const double diameter, + const double w, const double roughness, const double mu) + { + + double f = 0.; + // Reynolds number + const double re = std::abs(diameter * w / (area * mu)); + + if ( re == 0.0 ) { + // make sure it is because the mass rate is zero + assert(w == 0.); + return 0.0; + } + + const double re_value1 = 200.; + const double re_value2 = 4000.; + + if (re < re_value1) { + f = 16. / re; + } else if (re > re_value2){ + f = haalandFormular(re, diameter, roughness); + } else { // in between + const double f1 = 16. / re_value1; + const double f2 = haalandFormular(re_value2, diameter, roughness); + + f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1; + } + return f; + } + + + + + + + // calculating the friction pressure loss + // l is the segment length + // area is the segment cross area + // diameter is the segment inner diameter + // w is mass flow rate through the segment + // density is density + // roughness is the absolute roughness + // mu is the average phase viscosity + template + ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness, + const ValueType& density, const ValueType& w, const ValueType& mu) + { + const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value()); + // \Note: a factor of 2 needs to be here based on the dimensional analysis + return 2. * f * l * w * w / (area * area * diameter * density); + } + + + + + + template + ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density) + { + return (0.5 * mass_rate * mass_rate / (area * area * density)); + } + + +} // namespace mswellhelpers + +} + +#endif diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp new file mode 100644 index 000000000..fb44da239 --- /dev/null +++ b/opm/autodiff/MultisegmentWell.hpp @@ -0,0 +1,350 @@ +/* + 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_MULTISEGMENTWELL_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELL_HEADER_INCLUDED + + +#include + +namespace Opm +{ + + template + class MultisegmentWell: public WellInterface + { + public: + typedef WellInterface Base; + + using typename Base::WellState; + using typename Base::Simulator; + using typename Base::IntensiveQuantities; + using typename Base::FluidSystem; + using typename Base::ModelParameters; + using typename Base::MaterialLaw; + using typename Base::BlackoilIndices; + + /// the number of reservior equations + using Base::numEq; + + using Base::has_solvent; + using Base::has_polymer; + + // TODO: for now, not considering the polymer, solvent and so on to simplify the development process. + + // TODO: we need to have order for the primary variables and also the order for the well equations. + // sometimes, they are similar, while sometimes, they can have very different forms. + + // TODO: the following system looks not rather flexible. Looking into all kinds of possibilities + // TODO: gas is always there? how about oil water case? + // Is it gas oil two phase case? + static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0); + static const int GTotal = 0; + static const int WFrac = gasoil? -1000: 1; + static const int GFrac = gasoil? 1 : 2; + static const int SPres = gasoil? 2 : 3; + + /// the number of well equations // TODO: it should have a more general strategy for it + static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq : numEq + 1; + + using typename Base::Scalar; + using typename Base::ConvergenceReport; + + /// the matrix and vector types for the reservoir + 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; + + // TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell + // EvalR (Eval), EvalW, EvalRW + // TODO: for now, we only use one type to save some implementation efforts, while improve later. + typedef DenseAd::Evaluation EvalWell; + + MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param); + + 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& well_state) const; + + /// check whether the well equations get converged for this well + virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; + + /// 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, + 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); + + virtual void updatePrimaryVariables(const WellState& well_state) const; + + virtual void solveEqAndUpdateWellState(WellState& well_state); // const? + + virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& well_state); // should be const? + + /// number of segments for this well + /// int number_of_segments_; + int numberOfSegments() const; + + int numberOfPerforations() const; + + protected: + int number_segments_; + + // components of the pressure drop to be included + WellSegment::CompPressureDropEnum compPressureDrop() const; + // multi-phase flow model + WellSegment::MultiPhaseModelEnum multiphaseModel() const; + + // get the SegmentSet from the well_ecl_ + const SegmentSet& segmentSet() const; + + // protected member variables from the Base class + using Base::well_ecl_; + using Base::number_of_perforations_; // TODO: can use well_ecl_? + using Base::current_step_; + using Base::index_of_well_; + using Base::number_of_phases_; + + // TODO: the current implementation really relies on the order of the + // perforation does not change from the parser to Wells structure. + using Base::well_cells_; + using Base::param_; + using Base::well_index_; + using Base::well_type_; + using Base::first_perf_; + using Base::saturation_table_number_; + using Base::well_efficiency_factor_; + using Base::gravity_; + using Base::well_controls_; + using Base::perf_depth_; + + // protected functions from the Base class + using Base::active; + using Base::phaseUsage; + using Base::name; + using Base::numComponents; + using Base::flowPhaseToEbosPhaseIdx; + using Base::flowPhaseToEbosCompIdx; + using Base::getAllowCrossFlow; + + // TODO: trying to use the information from the Well opm-parser as much + // as possible, it will possibly be re-implemented later for efficiency reason. + + // the completions that is related to each segment + // the completions's ids are their index in the vector well_index_, well_cell_ + // This is also assuming the order of the completions in Well is the same with + // the order of the completions in wells. + // it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations + // belonging to this segment + std::vector > segment_perforations_; + + // the inlet segments for each segment. It is for convinience and efficiency reason + std::vector > segment_inlets_; + + // segment number is an ID of the segment, it is specified in the deck + // get the loation of the segment with a segment number in the segmentSet + int segmentNumberToIndex(const int segment_number) const; + + // TODO, the following should go to a class for computing purpose + // two off-diagonal matrices + mutable OffDiagMatWell duneB_; + mutable OffDiagMatWell duneC_; + // diagonal matrix for the well + mutable DiagMatWell duneD_; + + // residuals of the well equations + mutable BVectorWell resWell_; + + // 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_; + + // depth difference between perforations and the perforated grid cells + std::vector cell_perforation_depth_diffs_; + // pressure correction due to the different depth of the perforation and + // center depth of the grid block + std::vector cell_perforation_pressure_diffs_; + + // depth difference between the segment and the peforation + // or in another way, the depth difference between the perforation and + // the segment the perforation belongs to + std::vector perforation_segment_depth_diffs_; + + // the intial component compistion of segments + std::vector > segment_comp_initial_; + + // the densities of segment fluids + // we should not have this member variable + std::vector segment_densities_; + + // the viscosity of the segments + std::vector segment_viscosities_; + + // the mass rate of the segments + std::vector segment_mass_rates_; + + std::vector segment_depth_diffs_; + + void initMatrixAndVectors(const int num_cells) const; + + // protected functions + // EvalWell getBhp(); this one should be something similar to getSegmentPressure(); + // EvalWell getQs(); this one should be something similar to getSegmentRates() + // EValWell wellVolumeFractionScaled, wellVolumeFraction, wellSurfaceVolumeFraction ... these should have different names, and probably will be needed. + // bool crossFlowAllowed(const Simulator& ebosSimulator) const; probably will be needed + // 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 bool inner_iteration, + WellState& well_state) const; + + // initialize the segment rates with well rates + // when there is no more accurate way to initialize the segment rates, we initialize + // the segment rates based on well rates with a simple strategy + void initSegmentRatesWithWellRates(WellState& well_state) const; + + // computing the accumulation term for later use in well mass equations + void computeInitialComposition(); + + // compute the pressure difference between the perforation and cell center + void computePerfCellPressDiffs(const Simulator& ebosSimulator); + + // fraction value of the primary variables + // should we just use member variables to store them instead of calculating them again and again + EvalWell volumeFraction(const int seg, const int comp_idx) const; + + // F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p + EvalWell volumeFractionScaled(const int seg, const int comp_idx) const; + + // basically Q_p / \sigma_p Q_p + EvalWell surfaceVolumeFraction(const int seg, const int comp_idx) const; + + void computePerfRate(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s) const; + + // convert a Eval from reservoir to contain the derivative related to wells + EvalWell extendEval(const Eval& in) const; + + // compute the fluid properties, such as densities, viscosities, and so on, in the segments + // They will be treated implicitly, so they need to be of Evaluation type + void computeSegmentFluidProperties(const Simulator& ebosSimulator); + + EvalWell getSegmentPressure(const int seg) const; + + EvalWell getSegmentRate(const int seg, const int comp_idx) const; + + EvalWell getSegmentGTotal(const int seg) const; + + // get the mobility for specific perforation + void getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const; + + void assembleControlEq() const; + + void assemblePressureEq(const int seg) const; + + // hytrostatic pressure loss + EvalWell getHydroPressureLoss(const int seg) const; + + // frictinal pressure loss + EvalWell getFrictionPressureLoss(const int seg) const; + + void handleAccelerationPressureLoss(const int seg) const; + + // handling the overshooting and undershooting of the fractions + void processFractions(const int seg) const; + + void updateWellStateFromPrimaryVariables(WellState& well_state) const; + + double scalingFactor(const int comp_idx) const; + + bool frictionalPressureLossConsidered() const; + + bool accelerationalPressureLossConsidered() const; + + // TODO: try to make ebosSimulator const, as it should be + void iterateWellEquations(Simulator& ebosSimulator, + const double dt, + WellState& well_state); + + void assembleWellEqWithoutIteration(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells); + }; + +} + +#include "MultisegmentWell_impl.hpp" + +#endif // OPM_MULTISEGMENTWELL_HEADER_INCLUDED diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp new file mode 100644 index 000000000..c0310c0b2 --- /dev/null +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -0,0 +1,1893 @@ +/* + 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 + +namespace Opm +{ + + + template + MultisegmentWell:: + MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param) + : Base(well, time_step, wells, param) + , segment_perforations_(numberOfSegments()) + , segment_inlets_(numberOfSegments()) + , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) + , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0) + , perforation_segment_depth_diffs_(number_of_perforations_, 0.0) + , segment_comp_initial_(numberOfSegments(), std::vector(numComponents(), 0.0)) + , segment_densities_(numberOfSegments(), 0.0) + , segment_viscosities_(numberOfSegments(), 0.0) + , segment_mass_rates_(numberOfSegments(), 0.0) + , segment_depth_diffs_(numberOfSegments(), 0.0) + { + // not handling solvent or polymer for now with multisegment well + if (has_solvent) { + OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet"); + } + + if (has_polymer) { + OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet"); + } + // since we decide to use the SegmentSet from the well parser. we can reuse a lot from it. + // for other facilities needed but not available from parser, we need to process them here + + // initialize the segment_perforations_ + const CompletionSet& completion_set = well_ecl_->getCompletions(current_step_); + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const Completion& completion = completion_set.get(perf); + const int segment_number = completion.getSegmentNumber(); + const int segment_index = segmentNumberToIndex(segment_number); + segment_perforations_[segment_index].push_back(perf); + } + + // initialize the segment_inlets_ + for (int seg = 0; seg < numberOfSegments(); ++seg) { + const Segment& segment = segmentSet()[seg]; + const int segment_number = segment.segmentNumber(); + const int outlet_segment_number = segment.outletSegment(); + if (outlet_segment_number > 0) { + const int segment_index = segmentNumberToIndex(segment_number); + const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number); + segment_inlets_[outlet_segment_index].push_back(segment_index); + } + } + + // callcuate the depth difference between perforations and their segments + perf_depth_.resize(number_of_perforations_, 0.); + for (int seg = 0; seg < numberOfSegments(); ++seg) { + const double segment_depth = segmentSet()[seg].depth(); + for (const int perf : segment_perforations_[seg]) { + perf_depth_[perf] = completion_set.get(perf).getCenterDepth(); + perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth; + } + } + + // calculating the depth difference between the segment and its oulet_segments + // for the top segment, we will make its zero unless we find other purpose to use this value + for (int seg = 1; seg < numberOfSegments(); ++seg) { + const double segment_depth = segmentSet()[seg].depth(); + const int outlet_segment_number = segmentSet()[seg].outletSegment(); + const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)]; + const double outlet_depth = outlet_segment.depth(); + segment_depth_diffs_[seg] = segment_depth - outlet_depth; + } + } + + + + + + template + void + MultisegmentWell:: + 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); + + // TODO: for StandardWell, we need to update the perf depth here using depth_arg. + // for MultisegmentWell, it is much more complicated. + // It can be specified directly, it can be calculated from the segment depth, + // it can also use the cell center, which is the same for StandardWell. + // For the last case, should we update the depth with the depth_arg? For the + // future, it can be a source of wrong result with Multisegment well. + // An indicator from the opm-parser should indicate what kind of depth we should use here. + + // \Note: we do not update the depth here. And it looks like for now, we only have the option to use + // specified perforation depth + initMatrixAndVectors(num_cells); + + // calcuate the depth difference between the perforations and the perforated grid block + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf]; + } + } + + + + + + template + void + MultisegmentWell:: + initMatrixAndVectors(const int num_cells) const + { + duneB_.setBuildMode( OffDiagMatWell::row_wise ); + duneC_.setBuildMode( OffDiagMatWell::row_wise ); + duneD_.setBuildMode( DiagMatWell::row_wise ); + + // set the size and patterns for all the matrices and vectors + // [A C^T [x = [ res + // B D] x_well] res_well] + + // calculatiing the NNZ for duneD_ + // NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets) + { + int nnz_d = numberOfSegments(); + for (const std::vector& inlets : segment_inlets_) { + nnz_d += 2 * inlets.size(); + } + duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d); + } + duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_); + duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_); + + // we need to add the off diagonal ones + for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) { + // the number of the row corrspnds to the segment now + const int seg = row.index(); + // adding the item related to outlet relation + const Segment& segment = segmentSet()[seg]; + const int outlet_segment_number = segment.outletSegment(); + if (outlet_segment_number > 0) { // if there is a outlet_segment + const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number); + row.insert(outlet_segment_index); + } + + // Add nonzeros for diagonal + row.insert(seg); + + // insert the item related to its inlets + for (const int& inlet : segment_inlets_[seg]) { + row.insert(inlet); + } + } + + // make the C matrix + for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) { + // the number of the row corresponds to the segment number now. + for (const int& perf : segment_perforations_[row.index()]) { + const int cell_idx = well_cells_[perf]; + row.insert(cell_idx); + } + } + + // make the B^T matrix + for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) { + // the number of the row corresponds to the segment number now. + for (const int& perf : segment_perforations_[row.index()]) { + const int cell_idx = well_cells_[perf]; + row.insert(cell_idx); + } + } + + resWell_.resize( numberOfSegments() ); + + primary_variables_.resize(numberOfSegments()); + primary_variables_evaluation_.resize(numberOfSegments()); + } + + + + + + template + void + MultisegmentWell:: + initPrimaryVariablesEvaluation() const + { + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + primary_variables_evaluation_[seg][eq_idx] = 0.0; + primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]); + primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + numEq, 1.0); + } + } + } + + + + + + template + void + MultisegmentWell:: + assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) + { + + const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_; + if (use_inner_iterations) { + iterateWellEquations(ebosSimulator, dt, well_state); + } + + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells); + } + + + + + + template + void + MultisegmentWell:: + updateWellStateWithTarget(const int current, + WellState& well_state) const + { + // Updating well state bas on well control + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(well_controls_, current); + const double* distr = well_controls_iget_distr(well_controls_, current); + switch (well_controls_iget_type(well_controls_, current)) { + case BHP: { + well_state.bhp()[index_of_well_] = target; + const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_]; + // 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: { + well_state.thp()[index_of_well_] = target; + + /* const Opm::PhaseUsage& pu = phaseUsage(); + std::vector rates(3, 0.0); + if (active()[ Water ]) { + rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ]; + } + if (active()[ Oil ]) { + rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ]; + } + if (active()[ Gas ]) { + rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ]; + } */ + + // const int table_id = well_controls_iget_vfp(well_controls_, current); + // const double& thp = well_controls_iget_target(well_controls_, current); + // const double& alq = well_controls_iget_alq(well_controls_, current); + + // TODO: implement calculateBhpFromThp function + // well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current); + // also the top segment pressure + /* const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + well_state.segPress()[top_segment_index] = well_state.bhp()[index_of_well_]; */ + break; + } + + case RESERVOIR_RATE: // intentional fall-through + case SURFACE_RATE: + // for the update of the rates, after we update the well rates, we can try to scale + // the segment rates and perforation rates with the same factor + // or the other way, we can use the same approach like the initialization of the well state, + // we devide the well rates to update the perforation rates, then we update the segment rates based + // on the perforation rates. + // the second way is safer, since if the well control is changed, the first way will not guarantee the consistence + // of between the segment rates and peforation rates + + // checking the number of the phases under control + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < number_of_phases_; ++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 < number_of_phases_; ++phase) { + if (distr[phase] > 0.) { + well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target / distr[phase]; + } else { + well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = 0.; + } + } + + initSegmentRatesWithWellRates(well_state); + + } 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 < number_of_phases_; ++phase) { + if (distr[phase] > 0.0) { + original_rates_under_phase_control += well_state.wellRates()[number_of_phases_ * index_of_well_ + 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 < number_of_phases_; ++phase) { + well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] *= scaling_factor; + + // scaling the segment rates with the same way with well rates + const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + for (int seg = 0; seg < numberOfSegments(); ++seg) { + well_state.segRates()[number_of_phases_ * (seg + top_segment_index) + phase] *= scaling_factor; + } + } + } else { // scaling factor is not well defined 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 < number_of_phases_; ++phase) { + if (distr[phase] > 0.0) { + well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided / distr[phase]; + } else { + // this only happens for SURFACE_RATE control + well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] = target_rate_divided; + } + } + initSegmentRatesWithWellRates(well_state); + } + } + + break; + } // end of switch + + updatePrimaryVariables(well_state); + } + + + + + + template + void + MultisegmentWell:: + initSegmentRatesWithWellRates(WellState& well_state) const + { + for (int phase = 0; phase < number_of_phases_; ++phase) { + const double perf_phaserate = well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] / number_of_perforations_; + for (int perf = 0; perf < number_of_perforations_; ++perf) { + well_state.perfPhaseRates()[number_of_phases_ * (first_perf_ + perf) + phase] = perf_phaserate; + } + } + + const std::vector perforation_rates(well_state.perfPhaseRates().begin() + number_of_phases_ * first_perf_, + well_state.perfPhaseRates().begin() + + number_of_phases_ * (first_perf_ + number_of_perforations_) ); + std::vector segment_rates; + WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_, + 0, segment_rates); + const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + std::copy(segment_rates.begin(), segment_rates.end(), + well_state.segRates().begin() + number_of_phases_ * top_segment_index ); + // we need to check the top segment rates should be same with the well rates + } + + + + + + template + typename MultisegmentWell::ConvergenceReport + MultisegmentWell:: + getWellConvergence(const std::vector& B_avg) const + { + // assert((int(B_avg.size()) == numComponents()) || has_polymer); + assert( (int(B_avg.size()) == numComponents()) ); + + // checking if any residual is NaN or too large. The two large one is only handled for the well flux + std::vector> abs_residual(numberOfSegments(), std::vector(numWellEq, 0.0)); + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]); + } + } + + std::vector maximum_residual(numWellEq, 0.0); + + ConvergenceReport report; + // TODO: the following is a little complicated, maybe can be simplified in some way? + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + if (eq_idx < numComponents()) { // phase or component mass equations + const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx]; + // TODO: the report can not handle the segment number yet. + if (std::isnan(flux_residual)) { + report.nan_residual_found = true; + const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + report.nan_residual_wells.push_back(problem_well); + } else if (flux_residual > param_.max_residual_allowed_) { + report.too_large_residual_found = true; + const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + report.nan_residual_wells.push_back(problem_well); + } else { // it is a normal residual + if (flux_residual > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = flux_residual; + } + } + } else { // pressure equation + // TODO: we should distinguish the rate control equations, bhp control equations + // and the oridnary pressure equations + const double pressure_residal = abs_residual[seg][eq_idx]; + const std::string eq_name("Pressure"); + if (std::isnan(pressure_residal)) { + report.nan_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name}; + report.nan_residual_wells.push_back(problem_well); + } else if (std::isinf(pressure_residal)) { + report.too_large_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name}; + report.nan_residual_wells.push_back(problem_well); + } else { // it is a normal residual + if (pressure_residal > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = pressure_residal; + } + } + } + } + } + + if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found + // check convergence for flux residuals + for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) + { + report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_); + } + + report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_); + } else { // abnormal values found and no need to check the convergence + report.converged = false; + } + + return report; + } + + + + + + template + void + MultisegmentWell:: + apply(const BVector& x, BVector& Ax) const + { + BVectorWell Bx(duneB_.N()); + + duneB_.mv(x, Bx); + + // invDBx = duneD^-1 * Bx_ + const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx); + + // Ax = Ax - duneC_^T * invDBx + duneC_.mmtv(invDBx,Ax); + } + + + + + + template + void + MultisegmentWell:: + apply(BVector& r) const + { + // invDrw_ = duneD^-1 * resWell_ + const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_); + // r = r - duneC_^T * invDrw + duneC_.mmtv(invDrw, r); + } + + + + + + template + void + MultisegmentWell:: + recoverWellSolutionAndUpdateWellState(const BVector& x, + WellState& well_state) const + { + BVectorWell xw(1); + recoverSolutionWell(x, xw); + updateWellState(xw, false, well_state); + } + + + + + + template + void + MultisegmentWell:: + computeWellPotentials(const Simulator& /* ebosSimulator */, + const WellState& /* well_state */, + std::vector& /* well_potentials */) + { + OPM_THROW(std::runtime_error, "well potential calculation for multisegment wells is not supported yet"); + } + + + + + + template + void + MultisegmentWell:: + updatePrimaryVariables(const WellState& well_state) const + { + // TODO: to test using rate conversion coefficients to see if it will be better than + // this default one + + // the index of the top segment in the WellState + const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + const std::vector& segment_rates = well_state.segRates(); + const PhaseUsage& pu = phaseUsage(); + + for (int seg = 0; seg < numberOfSegments(); ++seg) { + // calculate the total rate for each segment + double total_seg_rate = 0.0; + const int seg_index = top_segment_index + seg; + // the segment pressure + primary_variables_[seg][SPres] = well_state.segPress()[seg_index]; + // TODO: under what kind of circustances, the following will be wrong? + // the definition of g makes the gas phase is always the last phase + for (int p = 0; p < number_of_phases_; p++) { + total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p]; + } + + primary_variables_[seg][GTotal] = total_seg_rate; + if (std::abs(total_seg_rate) > 0.) { + if (active()[Water]) { + const int water_pos = pu.phase_pos[Water]; + primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate; + } + if (active()[Gas]) { + const int gas_pos = pu.phase_pos[Gas]; + primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate; + } + } else { // total_seg_rate == 0 + if (well_type_ == INJECTOR) { + // only single phase injection handled + const double* distr = well_controls_get_current_distr(well_controls_); + if (active()[Water]) { + if (distr[pu.phase_pos[Water]] > 0.0) { + primary_variables_[seg][WFrac] = 1.0; + } else { + primary_variables_[seg][WFrac] = 0.0; + } + } + + if (active()[Gas]) { + if (distr[pu.phase_pos[Gas]] > 0.0) { + primary_variables_[seg][GFrac] = 1.0; + } else { + primary_variables_[seg][GFrac] = 0.0; + } + } + } else if (well_type_ == PRODUCER) { // producers + if (active()[Water]) { + primary_variables_[seg][WFrac] = 1.0 / number_of_phases_; + } + + if (active()[Gas]) { + primary_variables_[seg][GFrac] = 1.0 / number_of_phases_; + } + } + } + } + } + + + + + + template + void + MultisegmentWell:: + recoverSolutionWell(const BVector& x, BVectorWell& xw) const + { + BVectorWell resWell = resWell_; + // resWell = resWell - B * x + duneB_.mmv(x, resWell); + // xw = D^-1 * resWell + xw = mswellhelpers::invDXDirect(duneD_, resWell); + } + + + + + + template + void + MultisegmentWell:: + solveEqAndUpdateWellState(WellState& well_state) + { + // We assemble the well equations, then we check the convergence, + // which is why we do not put the assembleWellEq here. + const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + + updateWellState(dx_well, false, well_state); + } + + + + + + template + void + MultisegmentWell:: + computePerfCellPressDiffs(const Simulator& ebosSimulator) + { + for (int perf = 0; perf < number_of_perforations_; ++perf) { + + std::vector kr(number_of_phases_, 0.0); + std::vector density(number_of_phases_, 0.0); + + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + const auto& fs = intQuants.fluidState(); + + double sum_kr = 0.; + + const PhaseUsage& pu = phaseUsage(); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const int water_pos = pu.phase_pos[BlackoilPhases::Aqua]; + kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value(); + sum_kr += kr[water_pos]; + density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value(); + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oil_pos = pu.phase_pos[BlackoilPhases::Liquid]; + kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value(); + sum_kr += kr[oil_pos]; + density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value(); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + const int gas_pos = pu.phase_pos[BlackoilPhases::Vapour]; + kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value(); + sum_kr += kr[gas_pos]; + density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value(); + } + + assert(sum_kr != 0.); + + // calculate the average density + double average_density = 0.; + for (int p = 0; p < number_of_phases_; ++p) { + average_density += kr[p] * density[p]; + } + average_density /= sum_kr; + + cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf]; + } + } + + + + + + template + void + MultisegmentWell:: + computeInitialComposition() + { + for (int seg = 0; seg < numberOfSegments(); ++seg) { + // TODO: probably it should be numWellEq -1 more accurately, + // while by meaning it should be num_comp + for (int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) { + segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value(); + } + } + } + + + + + + template + void + MultisegmentWell:: + updateWellState(const BVectorWell& dwells, + const bool inner_iteration, + WellState& well_state) const + { + const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_; + + const double relaxation_factor = (use_inner_iterations && inner_iteration) ? 0.2 : 1.0; + + const double dFLimit = param_.dwell_fraction_max_; + const double max_pressure_change = param_.max_pressure_change_ms_wells_; + const std::vector > old_primary_variables = primary_variables_; + + for (int seg = 0; seg < numberOfSegments(); ++seg) { + if (active()[ Water ]) { + const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; + const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); + primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; + } + + if (active()[ Gas ]) { + const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; + const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); + primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; + } + + // handling the overshooting or undershooting of the fractions + processFractions(seg); + + // update the segment pressure + { + const int sign = dwells[seg][SPres] > 0.? 1 : -1; + const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); + primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited; + } + + // update the total rate // TODO: should we have a limitation of the total rate change? + { + primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal]; + } + + } + + updateWellStateFromPrimaryVariables(well_state); + } + + + + + + template + void + MultisegmentWell:: + calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& /* well_state */) + { + computePerfCellPressDiffs(ebosSimulator); + computeInitialComposition(); + } + + + + + + template + const SegmentSet& + MultisegmentWell:: + segmentSet() const + { + return well_ecl_->getSegmentSet(current_step_); + } + + + + + + template + int + MultisegmentWell:: + numberOfSegments() const + { + return segmentSet().numberSegment(); + } + + + + + + template + int + MultisegmentWell:: + numberOfPerforations() const + { + return segmentSet().number_of_perforations_; + } + + + + + + template + WellSegment::CompPressureDropEnum + MultisegmentWell:: + compPressureDrop() const + { + return segmentSet().compPressureDrop(); + } + + + + + + template + WellSegment::MultiPhaseModelEnum + MultisegmentWell:: + multiphaseModel() const + { + return segmentSet().multiPhaseModel(); + } + + + + + + template + int + MultisegmentWell:: + segmentNumberToIndex(const int segment_number) const + { + return segmentSet().segmentNumberToIndex(segment_number); + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + volumeFraction(const int seg, const int comp_idx) const + { + const PhaseUsage& pu = phaseUsage(); + + if (active()[Water] && comp_idx == pu.phase_pos[Water]) { + return primary_variables_evaluation_[seg][WFrac]; + } + + if (active()[Gas] && comp_idx == pu.phase_pos[Gas]) { + return primary_variables_evaluation_[seg][GFrac]; + } + + // Oil fraction + EvalWell oil_fraction = 1.0; + if (active()[Water]) { + oil_fraction -= primary_variables_evaluation_[seg][WFrac]; + } + + if (active()[Gas]) { + oil_fraction -= primary_variables_evaluation_[seg][GFrac]; + } + /* if (has_solvent) { + oil_fraction -= primary_variables_evaluation_[seg][SFrac]; + } */ + return oil_fraction; + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + volumeFractionScaled(const int seg, const int comp_idx) const + { + // For reservoir rate control, the distr in well control is used for the + // rate conversion coefficients. For the injection well, only the distr of the injection + // phase is not zero. + const double scale = scalingFactor(comp_idx); + if (scale > 0.) { + return volumeFraction(seg, comp_idx) / scale; + } + + return volumeFraction(seg, comp_idx); + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + surfaceVolumeFraction(const int seg, const int comp_idx) const + { + EvalWell sum_volume_fraction_scaled = 0.; + const int num_comp = numComponents(); + for (int idx = 0; idx < num_comp; ++idx) { + sum_volume_fraction_scaled += volumeFractionScaled(seg, idx); + } + + assert(sum_volume_fraction_scaled.value() != 0.); + + return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled; + } + + + + + + template + void + MultisegmentWell:: + computePerfRate(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s) const + { + const int num_comp = numComponents(); + std::vector cmix_s(num_comp, 0.0); + + // the composition of the components inside wellbore + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); + } + + const auto& fs = int_quants.fluidState(); + + const EvalWell pressure_cell = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + const EvalWell rs = extendEval(fs.Rs()); + const EvalWell rv = extendEval(fs.Rv()); + + // not using number_of_phases_ because of solvent + std::vector b_perfcells(num_comp, 0.0); + + for (int phase = 0; phase < number_of_phases_; ++phase) { + const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase); + b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos)); + } + + // pressure difference between the segment and the perforation + const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf]; + // pressure difference between the perforation and the grid cell + const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf]; + + // Pressure drawdown (also used to determine direction of flow) + // TODO: not 100% sure about the sign of the seg_perf_press_diff + const EvalWell drawdown = (pressure_cell + cell_perf_press_diff) - (segment_pressure + perf_seg_press_diff); + + const Opm::PhaseUsage& pu = phaseUsage(); + + // producing perforations + if ( drawdown > 0.0) { + // Do nothing is crossflow is not allowed + if (!allow_cf && well_type_ == INJECTOR) { + return; + } + + // compute component volumetric rates at standard conditions + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown); + cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; + } + + if (active()[Oil] && active()[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const EvalWell cq_s_oil = cq_s[oilpos]; + const EvalWell cq_s_gas = cq_s[gaspos]; + cq_s[gaspos] += rs * cq_s_oil; + cq_s[oilpos] += rv * cq_s_gas; + } + } else { // injecting perforations + // Do nothing if crossflow is not allowed + if (!allow_cf && well_type_ == PRODUCER) { + return; + } + + // for injecting perforations, we use total mobility + EvalWell total_mob = mob_perfcells[0]; + for (int comp_idx = 1; comp_idx < num_comp; ++comp_idx) { + total_mob += mob_perfcells[comp_idx]; + } + + // injection perforations total volume rates + const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown); + + // compute volume ratio between connection and at standard conditions + EvalWell volume_ratio = 0.0; + if (active()[Water]) { + const int watpos = pu.phase_pos[Water]; + volume_ratio += cmix_s[watpos] / b_perfcells[watpos]; + } + + 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 + // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations + // basically, for injecting perforations, the wellbore is the upstreaming side. + 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; + volume_ratio += tmp_oil / b_perfcells[oilpos]; + + const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d; + volume_ratio += tmp_gas / b_perfcells[gaspos]; + } else { // not having gas and oil at the same time + if (active()[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + volume_ratio += cmix_s[oilpos] / b_perfcells[oilpos]; + } + if (active()[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + volume_ratio += cmix_s[gaspos] / b_perfcells[gaspos]; + } + } + // injecting connections total volumerates at standard conditions + EvalWell cqt_is = cqt_i / volume_ratio; + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; + } + } // end for injection perforations + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + extendEval(const Eval& in) const + { + EvalWell out = 0.0; + out.setValue(in.value()); + for(int eq_idx = 0; eq_idx < numEq;++eq_idx) { + out.setDerivative(eq_idx, in.derivative(eq_idx)); + } + return out; + } + + + + + + template + void + MultisegmentWell:: + computeSegmentFluidProperties(const Simulator& ebosSimulator) + { + // TODO: the concept of phases and components are rather confusing in this function. + // needs to be addressed sooner or later. + + // get the temperature for later use. It is only useful when we are not handling + // thermal related simulation + // basically, it is a single value for all the segments + + EvalWell temperature; + // not sure how to handle the pvt region related to segment + // for the current approach, we use the pvt region of the first perforated cell + // although there are some text indicating using the pvt region of the lowest + // perforated cell + // TODO: later to investigate how to handle the pvt region + int pvt_region_index; + { + // using the first perforated cell + const int cell_idx = well_cells_[0]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value()); + pvt_region_index = fs.pvtRegionIndex(); + } + + std::vector surf_dens(number_of_phases_); + // Surface density. + // not using num_comp here is because solvent can be component + for (int phase = 0; phase < number_of_phases_; ++phase) { + surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index ); + } + + const int num_comp = numComponents(); + const Opm::PhaseUsage& pu = phaseUsage(); + for (int seg = 0; seg < numberOfSegments(); ++seg) { + // the compostion of the components inside wellbore under surface condition + std::vector mix_s(num_comp, 0.0); + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); + } + + std::vector b(num_comp, 0.0); + // it is the phase viscosities asked for + std::vector visc(number_of_phases_, 0.0); + const EvalWell seg_pressure = getSegmentPressure(seg); + if (pu.phase_used[BlackoilPhases::Aqua]) { + // TODO: what is the difference between Water and BlackoilPhases::Aqua? + const int water_pos = pu.phase_pos[BlackoilPhases::Aqua]; + b[water_pos] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[water_pos] = + FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure); + } + + EvalWell rv(0.0); + // gas phase + if (pu.phase_used[BlackoilPhases::Vapour]) { + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour]; + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid]; + const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[oilpos] > 0.0) { + if (mix_s[gaspos] > 0.0) { + rv = mix_s[oilpos] / mix_s[gaspos]; + } + + if (rv > rvmax) { + rv = rvmax; + } + b[gaspos] = + FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); + visc[gaspos] = + FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv); + } else { // no oil exists + b[gaspos] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[gaspos] = + FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[gaspos] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[gaspos] = + FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } + + EvalWell rs(0.0); + // oil phase + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid]; + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour]; + const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[gaspos] > 0.0) { + if (mix_s[oilpos] > 0.0) { + rs = mix_s[gaspos] / mix_s[oilpos]; + } + + if (rs > rsmax) { + rs = rsmax; + } + b[oilpos] = + FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); + visc[oilpos] = + FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs); + } else { // no oil exists + b[oilpos] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[oilpos] = + FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[oilpos] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[oilpos] = + FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } + + std::vector mix(mix_s); + if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour]; + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid]; + if (rs != 0.0) { // rs > 0.0? + mix[gaspos] = (mix_s[gaspos] - mix_s[oilpos] * rs) / (1. - rs * rv); + } + if (rv != 0.0) { // rv > 0.0? + mix[oilpos] = (mix_s[oilpos] - mix_s[gaspos] * rv) / (1. - rs * rv); + } + } + + EvalWell volrat(0.0); + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + volrat += mix[comp_idx] / b[comp_idx]; + } + + segment_viscosities_[seg] = 0.; + // calculate the average viscosity + for (int p = 0; p < number_of_phases_; ++p) { + const EvalWell phase_fraction = mix[p] / b[p] / volrat; + segment_viscosities_[seg] += visc[p] * phase_fraction; + } + + EvalWell density(0.0); + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + density += surf_dens[comp_idx] * mix_s[comp_idx]; + } + segment_densities_[seg] = density / volrat; + + // calculate the mass rates + segment_mass_rates_[seg] = 0.; + for (int phase = 0; phase < number_of_phases_; ++phase) { + const EvalWell rate = getSegmentRate(seg, phase); + segment_mass_rates_[seg] += rate * surf_dens[phase]; + } + } + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getSegmentPressure(const int seg) const + { + return primary_variables_evaluation_[seg][SPres]; + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getSegmentRate(const int seg, + const int comp_idx) const + { + return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx); + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getSegmentGTotal(const int seg) const + { + return primary_variables_evaluation_[seg][GTotal]; + } + + + + + + template + void + MultisegmentWell:: + getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const + { + // TODO: most of this function, if not the whole function, can be moved to the base class + 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)); + } + + } + } + + + + + + template + void + MultisegmentWell:: + assembleControlEq() const + { + EvalWell control_eq(0.0); + + switch (well_controls_get_current_type(well_controls_)) { + case THP: // not handling this one for now + { + OPM_THROW(std::runtime_error, "Not handling THP control for Multisegment wells for now"); + } + case BHP: + { + const double target_bhp = well_controls_get_current_target(well_controls_); + control_eq = getSegmentPressure(0) - target_bhp; + break; + } + case SURFACE_RATE: + { + // finding if it is a single phase control or combined phase control + int number_phases_under_control = 0; + const double* distr = well_controls_get_current_distr(well_controls_); + for (int phase = 0; phase < number_of_phases_; ++phase) { + if (distr[phase] > 0.0) { + ++number_phases_under_control; + } + } + assert(number_phases_under_control > 0); + + const std::vector g = {1.0, 1.0, 0.01}; + const double target_rate = well_controls_get_current_target(well_controls_); + // TODO: the two situations below should be able to be merged to be handled as one situation + + if (number_phases_under_control == 1) { // single phase control + for (int phase = 0; phase < number_of_phases_; ++phase) { + if (distr[phase] > 0.) { // under the control of this phase + control_eq = getSegmentGTotal(0) * volumeFraction(0, phase) - g[phase] * target_rate; + break; + } + } + } else { // multiphase rate control + EvalWell rate_for_control(0.0); + const EvalWell G_total = getSegmentGTotal(0); + for (int phase = 0; phase < number_of_phases_; ++phase) { + if (distr[phase] > 0.) { + rate_for_control += G_total * volumeFractionScaled(0, phase); + } + } + // TODO: maybe the following equation can be scaled a little bit for gas phase + control_eq = rate_for_control - target_rate; + } + break; + } + case RESERVOIR_RATE: + { + EvalWell rate_for_control(0.0); // reservoir rate + const double* distr = well_controls_get_current_distr(well_controls_); + for (int phase = 0; phase < number_of_phases_; ++phase) { + if (distr[phase] > 0.) { + rate_for_control += getSegmentGTotal(0) * volumeFraction(0, phase); + } + } + const double target_rate = well_controls_get_current_target(well_controls_); + control_eq = rate_for_control - target_rate; + break; + } + default: + OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name()); + } + + + // using control_eq to update the matrix and residuals + + resWell_[0][SPres] = control_eq.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq); + } + } + + + + + + template + void + MultisegmentWell:: + assemblePressureEq(const int seg) const + { + assert(seg != 0); // not top segment + + // for top segment, the well control equation will be used. + EvalWell pressure_equation = getSegmentPressure(seg); + + // we need to handle the pressure difference between the two segments + // we only consider the hydrostatic pressure loss first + pressure_equation -= getHydroPressureLoss(seg); + + if (frictionalPressureLossConsidered()) { + pressure_equation -= getFrictionPressureLoss(seg); + } + + resWell_[seg][SPres] = pressure_equation.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][seg][SPres][pv_idx] = pressure_equation.derivative(pv_idx + numEq); + } + + // contribution from the outlet segment + const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); + + resWell_[seg][SPres] -= outlet_pressure.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); + } + + if (accelerationalPressureLossConsidered()) { + handleAccelerationPressureLoss(seg); + } + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getHydroPressureLoss(const int seg) const + { + return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg]; + } + + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getFrictionPressureLoss(const int seg) const + { + const EvalWell mass_rate = segment_mass_rates_[seg]; + const EvalWell density = segment_densities_[seg]; + const EvalWell visc = segment_viscosities_[seg]; + const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength(); + assert(length > 0.); + const double roughness = segmentSet()[seg].roughness(); + const double area = segmentSet()[seg].crossArea(); + const double diameter = segmentSet()[seg].internalDiameter(); + + const double sign = mass_rate < 0. ? 1.0 : - 1.0; + + return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc); + } + + + + + + template + void + MultisegmentWell:: + handleAccelerationPressureLoss(const int seg) const + { + // TODO: this pressure loss is not significant enough to be well tested yet. + // handle the out velcocity head + const double area = segmentSet()[seg].crossArea(); + const EvalWell mass_rate = segment_mass_rates_[seg]; + const EvalWell density = segment_densities_[seg]; + const EvalWell out_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density); + + resWell_[seg][SPres] -= out_velocity_head.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][seg][SPres][pv_idx] -= out_velocity_head.derivative(pv_idx + numEq); + } + + // calcuate the maximum cross-area among the segment and its inlet segments + double max_area = area; + for (const int inlet : segment_inlets_[seg]) { + const double inlet_area = segmentSet()[inlet].crossArea(); + if (inlet_area > max_area) { + max_area = inlet_area; + } + } + + // handling the velocity head of intlet segments + for (const int inlet : segment_inlets_[seg]) { + const EvalWell density = segment_densities_[inlet]; + const EvalWell mass_rate = segment_mass_rates_[inlet]; + const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density); + resWell_[seg][SPres] += inlet_velocity_head.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq); + } + } + } + + + + + + template + void + MultisegmentWell:: + processFractions(const int seg) const + { + const PhaseUsage& pu = phaseUsage(); + + std::vector fractions(number_of_phases_, 0.0); + + assert( active()[Oil] ); + const int oil_pos = pu.phase_pos[Oil]; + fractions[oil_pos] = 1.0; + + if ( active()[Water] ) { + const int water_pos = pu.phase_pos[Water]; + fractions[water_pos] = primary_variables_[seg][WFrac]; + fractions[oil_pos] -= fractions[water_pos]; + } + + if ( active()[Gas] ) { + const int gas_pos = pu.phase_pos[Gas]; + fractions[gas_pos] = primary_variables_[seg][GFrac]; + fractions[oil_pos] -= fractions[gas_pos]; + } + + if (active()[Water]) { + const int water_pos = pu.phase_pos[Water]; + if (fractions[water_pos] < 0.0) { + if ( active()[Gas] ) { + fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]); + } + fractions[oil_pos] /= (1.0 - fractions[water_pos]); + fractions[water_pos] = 0.0; + } + } + + if (active()[Gas]) { + const int gas_pos = pu.phase_pos[Gas]; + if (fractions[gas_pos] < 0.0) { + if ( active()[Water] ) { + fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]); + } + fractions[oil_pos] /= (1.0 - fractions[gas_pos]); + fractions[gas_pos] = 0.0; + } + } + + if (fractions[oil_pos] < 0.0) { + if ( active()[Water] ) { + fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]); + } + if ( active()[Gas] ) { + fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]); + } + fractions[oil_pos] = 0.0; + } + + if ( active()[Water] ) { + primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]]; + } + + if ( active()[Gas] ) { + primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]]; + } + } + + + + + + template + void + MultisegmentWell:: + updateWellStateFromPrimaryVariables(WellState& well_state) const + { + const PhaseUsage& pu = phaseUsage(); + assert( active()[Oil] ); + const int oil_pos = pu.phase_pos[Oil]; + + for (int seg = 0; seg < numberOfSegments(); ++seg) { + std::vector fractions(number_of_phases_, 0.0); + fractions[oil_pos] = 1.0; + + if ( active()[Water] ) { + const int water_pos = pu.phase_pos[Water]; + fractions[water_pos] = primary_variables_[seg][WFrac]; + fractions[oil_pos] -= fractions[water_pos]; + } + + if ( active()[Gas] ) { + const int gas_pos = pu.phase_pos[Gas]; + fractions[gas_pos] = primary_variables_[seg][GFrac]; + fractions[oil_pos] -= fractions[gas_pos]; + } + + // convert the fractions to be Q_p / G_total to calculate the phase rates + for (int p = 0; p < number_of_phases_; ++p) { + const double scale = scalingFactor(p); + // for injection wells, there should only one non-zero scaling factor + if (scale > 0.) { + fractions[p] /= scale; + } else { + // this should only happens to injection wells + fractions[p] = 0.; + } + } + + // calculate the phase rates based on the primary variables + const double g_total = primary_variables_[seg][GTotal]; + const int top_segment_index = well_state.topSegmentIndex(index_of_well_); + for (int p = 0; p < number_of_phases_; ++p) { + const double phase_rate = g_total * fractions[p]; + well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate; + if (seg == 0) { // top segment + well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate; + } + } + + // update the segment pressure + well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres]; + if (seg == 0) { // top segment + well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index]; + } + } + } + + + + + + template + double + MultisegmentWell:: + scalingFactor(const int comp_idx) const + { + const double* distr = well_controls_get_current_distr(well_controls_); + + if (well_controls_get_current_type(well_controls_) == RESERVOIR_RATE) { + // if (has_solvent && phaseIdx == contiSolventEqIdx ) + // OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented"); + return distr[comp_idx]; + } + + const PhaseUsage& pu = phaseUsage(); + + if (active()[Water] && pu.phase_pos[Water] == comp_idx) + return 1.0; + if (active()[Oil] && pu.phase_pos[Oil] == comp_idx) + return 1.0; + if (active()[Gas] && pu.phase_pos[Gas] == comp_idx) + return 0.01; + // if (has_solvent && phaseIdx == contiSolventEqIdx ) + // return 0.01; + + // we should not come this far + assert(false); + return 1.0; + } + + + + + + template + bool + MultisegmentWell:: + frictionalPressureLossConsidered() const + { + // HF- and HFA needs to consider frictional pressure loss + return (segmentSet().compPressureDrop() != WellSegment::H__); + } + + + + + + template + bool + MultisegmentWell:: + accelerationalPressureLossConsidered() const + { + return (segmentSet().compPressureDrop() == WellSegment::HFA); + } + + + + + + template + void + MultisegmentWell:: + iterateWellEquations(Simulator& ebosSimulator, + const double dt, + WellState& well_state) + { + // basically, it only iterate through the equations. + // we update the primary variables + // if converged, we can update the well_state. + // the function updateWellState() should have a flag to show + // if we will update the well state. + const int max_iter_number = param_.max_inner_iter_ms_wells_; + int it = 0; + for (; it < max_iter_number; ++it) { + + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true); + + const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + + // TODO: use these small values for now, not intend to reach the convergence + // in this stage, but, should we? + // We should try to avoid hard-code values in the code. + // If we want to use the real one, we need to find a way to get them. + // const std::vector B {0.8, 0.8, 0.008}; + const std::vector B {0.5, 0.5, 0.005}; + + const ConvergenceReport report = getWellConvergence(B); + + if (report.converged) { + break; + } + + updateWellState(dx_well, true, well_state); + + initPrimaryVariablesEvaluation(); + } + // TODO: maybe we should not use these values if they are not converged. + } + + + + + + template + void + MultisegmentWell:: + assembleWellEqWithoutIteration(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) + { + // calculate the fluid properties needed. + computeSegmentFluidProperties(ebosSimulator); + + // clear all entries + if (!only_wells) { + duneB_ = 0.0; + duneC_ = 0.0; + } + + duneD_ = 0.0; + resWell_ = 0.0; + + // for the black oil cases, there will be four equations, + // the first three of them are the mass balance equations, the last one is the pressure equations. + // + // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same. + + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + const bool allow_cf = getAllowCrossFlow(); + + const int nseg = numberOfSegments(); + const int num_comp = numComponents(); + + for (int seg = 0; seg < nseg; ++seg) { + // calculating the accumulation term // TODO: without considering the efficiencty factor for now + // volume of the segment + { + const double volume = segmentSet()[seg].volume(); + // for each component + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx]) + + getSegmentRate(seg, comp_idx); + + resWell_[seg][comp_idx] += accumulation_term.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq); + } + } + } + + // considering the contributions from the inlet segments + { + for (const int inlet : segment_inlets_[seg]) { + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx); + resWell_[seg][comp_idx] -= inlet_rate.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq); + } + } + } + } + + // calculating the perforation rate for each perforation that belongs to this segment + const EvalWell seg_pressure = getSegmentPressure(seg); + for (const int perf : segment_perforations_[seg]) { + const int cell_idx = well_cells_[perf]; + const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector mob(num_comp, 0.0); + getMobility(ebosSimulator, perf, mob); + std::vector cq_s(num_comp, 0.0); + computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s); + + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + // the cq_s entering mass balance equations need to consider the efficiency factors. + const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_; + + if (!only_wells) { + // subtract sum of component fluxes in the reservoir equation. + // need to consider the efficiency factor + // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input + // is a component index :D + ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value(); + } + + // subtract sum of phase fluxes in the well equations. + resWell_[seg][comp_idx] -= cq_s_effective.value(); + + // assemble the jacobians + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix + } + // the index name for the D should be eq_idx / pv_idx + duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq); + } + + for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx); + duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); + } + } + } + // TODO: we should save the perforation pressure and preforation rates? + // we do not use it in the simulation for now, while we might need them if + // we handle the pressure in SEG mode. + } + + // the fourth dequation, the pressure drop equation + if (seg == 0) { // top segment, pressure equation is the control equation + assembleControlEq(); + } else { + assemblePressureEq(seg); + } + } + } + +} diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index adda1ebf7..f6fadeeb7 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -66,7 +66,7 @@ public: typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; - typedef StandardWellsDense WellModel; + typedef BlackoilWellModel WellModel; typedef RateConverter::SurfaceToReservoirVoidage > RateConverterType; @@ -301,6 +301,20 @@ public: WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, rateConverter_, terminal_output_, timer.currentStepNum(), legacyCellPvtRegionIdx_); + // handling MS well related + if (model_param_.use_multisegment_well_) { // if we use MultisegmentWell model + for (const auto& well : wells_ecl) { + // TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well + // while keep the well shut. More accurately, we should check if the well exisits in the Wells + // structure here + if (well->isMultiSegment(timer.currentStepNum()) ) { // there is one well is MS well + well_state.initWellStateMSWell(wells, wells_ecl, timer.currentStepNum(), phaseUsage_, prev_well_state); + break; + } + } + } + + auto solver = createSolver(well_model); std::vector> currentFluidInPlace; @@ -766,9 +780,9 @@ protected: return totals; } - + void outputTimestampFIP(SimulatorTimer& timer, const std::string version) - { + { std::ostringstream ss; boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d %b %Y"); ss.imbue(std::locale(std::locale::classic(), facet)); diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index b32b25344..eec5c8396 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -99,7 +99,7 @@ namespace Opm static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; - StandardWell(const Well* well, const int time_step, const Wells* wells); + StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector* active_arg, @@ -120,21 +120,8 @@ namespace Opm 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); + virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; /// Ax = Ax - C D^-1 B x virtual void apply(const BVector& x, BVector& Ax) const; @@ -143,16 +130,20 @@ namespace Opm /// 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, + virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, 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; + std::vector& well_potentials) /* const */; virtual void updatePrimaryVariables(const WellState& well_state) const; + virtual void solveEqAndUpdateWellState(WellState& well_state); + + virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& well_state); // should be const? protected: // protected functions from the Base class @@ -170,8 +161,8 @@ namespace Opm // protected member variables from the Base class using Base::vfp_properties_; using Base::gravity_; + using Base::param_; using Base::well_efficiency_factor_; - using Base::phase_usage_; using Base::first_perf_; using Base::ref_depth_; using Base::perf_depth_; @@ -240,7 +231,6 @@ namespace Opm // 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 @@ -268,8 +258,11 @@ namespace Opm const std::vector& rvmax_perf, const std::vector& surf_dens_perf); - virtual void solveEqAndUpdateWellState(const ModelParameters& param, - WellState& well_state); + // computing the accumulation term for later use in well mass equations + void computeAccumWell(); + + void computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw); // TODO: to check whether all the paramters are required void computePerfRate(const IntensiveQuantities& intQuants, diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index a596ccb75..6bd6bd1bd 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -24,8 +24,8 @@ namespace Opm { template StandardWell:: - StandardWell(const Well* well, const int time_step, const Wells* wells) - : Base(well, time_step, wells) + StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param) + : Base(well, time_step, wells, param) , perf_densities_(number_of_perforations_) , perf_pressure_diffs_(number_of_perforations_) , primary_variables_(numWellEq, 0.0) @@ -73,7 +73,6 @@ namespace Opm } 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); @@ -317,7 +316,7 @@ namespace Opm StandardWell:: wellVolumeFraction(const int compIdx) const { - const auto pu = phaseUsage(); + const auto& pu = phaseUsage(); if (active()[Water] && compIdx == pu.phase_pos[Water]) { return primary_variables_evaluation_[WFrac]; } @@ -401,14 +400,14 @@ namespace Opm for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); } - auto& fs = intQuants.fluidState(); + const auto& fs = intQuants.fluidState(); - EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); - EvalWell rs = extendEval(fs.Rs()); - EvalWell rv = extendEval(fs.Rv()); + const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + const EvalWell rs = extendEval(fs.Rs()); + const EvalWell rv = extendEval(fs.Rv()); std::vector b_perfcells_dense(numComp, 0.0); for (int phase = 0; phase < np; ++phase) { - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); } if (has_solvent) { @@ -416,8 +415,8 @@ namespace Opm } // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + cdp; - EvalWell drawdown = pressure - well_pressure; + const EvalWell well_pressure = bhp + cdp; + const EvalWell drawdown = pressure - well_pressure; // producing perforations if ( drawdown.value() > 0 ) { @@ -523,8 +522,10 @@ namespace Opm const int np = number_of_phases_; // clear all entries - duneB_ = 0.0; - duneC_ = 0.0; + if (!only_wells) { + duneB_ = 0.0; + duneC_ = 0.0; + } invDuneD_ = 0.0; resWell_ = 0.0; @@ -763,12 +764,11 @@ namespace Opm void StandardWell:: updateWellState(const BVectorWell& dwells, - const BlackoilModelParameters& param, WellState& well_state) const { const int np = number_of_phases_; - const double dBHPLimit = param.dbhp_max_rel_; - const double dFLimit = param.dwell_fraction_max_; + const double dBHPLimit = param_.dbhp_max_rel_; + const double dFLimit = param_.dwell_fraction_max_; const auto pu = phaseUsage(); const std::vector xvar_well_old = primary_variables_; @@ -1114,69 +1114,6 @@ namespace Opm - template - void - StandardWell:: - updateWellControl(WellState& xw, - wellhelpers::WellSwitchingLogger& logger) const - { - const int np = number_of_phases_; - 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:: @@ -1190,7 +1127,7 @@ namespace Opm const int nperf = number_of_perforations_; // TODO: can make this a member? const int numComp = numComponents(); - const PhaseUsage& pu = *phase_usage_; + const PhaseUsage& pu = phaseUsage(); b_perf.resize(nperf*numComp); surf_dens_perf.resize(nperf*numComp); const int w = index_of_well_; @@ -1298,7 +1235,7 @@ namespace Opm const int np = number_of_phases_; const int nperf = number_of_perforations_; const int num_comp = numComponents(); - const PhaseUsage* phase_usage = phase_usage_; + const PhaseUsage& phase_usage = phaseUsage(); // 1. Compute the flow (in surface volume units for each // component) exiting up the wellbore from each perforation, @@ -1326,8 +1263,8 @@ namespace Opm // 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]; + 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); @@ -1426,9 +1363,7 @@ namespace Opm template typename StandardWell::ConvergenceReport StandardWell:: - getWellConvergence(Simulator& ebosSimulator, - const std::vector& B_avg, - const ModelParameters& param) const + getWellConvergence(const std::vector& B_avg) const { typedef double Scalar; typedef std::vector< Scalar > Vector; @@ -1440,8 +1375,8 @@ namespace Opm // 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_; + 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. @@ -1553,15 +1488,28 @@ namespace Opm template void StandardWell:: - solveEqAndUpdateWellState(const ModelParameters& param, - WellState& well_state) + solveEqAndUpdateWellState(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); + updateWellState(dx_well, well_state); + } + + + + + + template + void + StandardWell:: + calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& well_state) + { + computeWellConnectionPressures(ebosSimulator, well_state); + computeAccumWell(); } @@ -1573,6 +1521,8 @@ namespace Opm StandardWell:: computeAccumWell() { + // TODO: it should be num_comp, while it also bring problem for + // the polymer case. for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value(); } @@ -1642,12 +1592,11 @@ namespace Opm void StandardWell:: recoverWellSolutionAndUpdateWellState(const BVector& x, - const ModelParameters& param, WellState& well_state) const { BVectorWell xw(1); recoverSolutionWell(x, xw); - updateWellState(xw, param, well_state); + updateWellState(xw, well_state); } @@ -1723,7 +1672,7 @@ namespace Opm 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_; + const Opm::PhaseUsage& pu = phaseUsage(); std::vector rates(3, 0.0); if (active()[ Water ]) { @@ -1794,10 +1743,16 @@ namespace Opm StandardWell:: computeWellPotentials(const Simulator& ebosSimulator, const WellState& well_state, - std::vector& well_potentials) const + std::vector& well_potentials) // const { - const int np = number_of_phases_; + 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(); + + const int np = number_of_phases_; well_potentials.resize(np, 0.0); // get the bhp value based on the bhp constraints diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index c78875214..8f33f54fb 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -86,7 +86,7 @@ namespace Opm static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); /// Constructor - WellInterface(const Well* well, const int time_step, const Wells* wells); + WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param); /// Virutal destructor virtual ~WellInterface() {} @@ -144,12 +144,9 @@ namespace Opm } }; - virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, - const std::vector& B_avg, - const ModelParameters& param) const = 0; + virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const = 0; - virtual void solveEqAndUpdateWellState(const ModelParameters& param, - WellState& well_state) = 0; + virtual void solveEqAndUpdateWellState(WellState& well_state) = 0; virtual void assembleWellEq(Simulator& ebosSimulator, const double dt, @@ -165,7 +162,7 @@ namespace Opm /// 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, + virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const = 0; /// Ax = Ax - C D^-1 B x @@ -174,25 +171,22 @@ namespace Opm /// r = r - C D^-1 Rw virtual void apply(BVector& r) const = 0; + // TODO: before we decide to put more information under mutable, this function is not const 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; + std::vector& well_potentials) = 0; virtual void updateWellStateWithTarget(const int current, WellState& xw) const = 0; - virtual void updateWellControl(WellState& xw, - wellhelpers::WellSwitchingLogger& logger) const = 0; + void updateWellControl(WellState& xw, + wellhelpers::WellSwitchingLogger& logger) const; virtual void updatePrimaryVariables(const WellState& well_state) const = 0; + virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, + const WellState& xw) = 0; // should be const? + protected: // to indicate a invalid connection @@ -205,6 +199,9 @@ namespace Opm // the index of well in Wells struct int index_of_well_; + // simulation parameters + const ModelParameters& param_; + // well type // INJECTOR or PRODUCER enum WellType well_type_; @@ -217,21 +214,18 @@ namespace Opm 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_; @@ -273,7 +267,7 @@ namespace Opm int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - // TODO: it is dumplicated with StandardWellsDense + // TODO: it is dumplicated with BlackoilWellModel int numComponents() const; double wsolvent() const; diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 8e5107bc5..7f5122070 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -25,9 +25,10 @@ namespace Opm template WellInterface:: - WellInterface(const Well* well, const int time_step, const Wells* wells) + WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param) : well_ecl_(well) , current_step_(time_step) + , param_(param) { if (!well) { OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface"); @@ -401,13 +402,76 @@ namespace Opm + template + void + WellInterface:: + updateWellControl(WellState& well_state, + wellhelpers::WellSwitchingLogger& logger) const + { + const int np = number_of_phases_; + const int w = index_of_well_; + + const int old_control_index = well_state.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 = well_state.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( + well_state.bhp(), well_state.thp(), well_state.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. + well_state.currentControls()[w] = ctrl_index; + current = well_state.currentControls()[w]; + well_controls_set_current( wc, current); + } + + // the new well control indices after all the related updates, + const int updated_control_index = well_state.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, well_state); + } + } + + + + + template bool WellInterface:: checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, const WellState& well_state) const { - const Opm::PhaseUsage& pu = *phase_usage_; + const Opm::PhaseUsage& pu = phaseUsage(); const int np = number_of_phases_; if (econ_production_limits.onMinOilRate()) { @@ -464,7 +528,7 @@ namespace Opm double violation_extent = -1.0; const int np = number_of_phases_; - const Opm::PhaseUsage& pu = *phase_usage_; + const Opm::PhaseUsage& pu = phaseUsage(); const int well_number = index_of_well_; assert(active()[Oil]); diff --git a/opm/autodiff/WellMultiSegment.cpp b/opm/autodiff/WellMultiSegment.cpp index 5e8e4704c..2d970040f 100644 --- a/opm/autodiff/WellMultiSegment.cpp +++ b/opm/autodiff/WellMultiSegment.cpp @@ -60,7 +60,7 @@ namespace Opm // we change the ID to location now for easier use later. for (int i = 0; i < m_number_of_segments_; ++i) { // The segment number for top segment is 0, the segment number of its outlet segment will be -1 - m_outlet_segment_[i] = segment_set.numberToLocation(segment_set[i].outletSegment()); + m_outlet_segment_[i] = segment_set.segmentNumberToIndex(segment_set[i].outletSegment()); m_segment_length_[i] = segment_set[i].totalLength(); m_segment_depth_[i] = segment_set[i].depth(); m_segment_internal_diameter_[i] = segment_set[i].internalDiameter(); @@ -111,7 +111,7 @@ namespace Opm int i_segment = completion_set.get(i).getSegmentNumber(); // using the location of the segment in the array as the segment number/id. // TODO: it can be helpful for output or postprocessing if we can keep the original number. - i_segment = segment_set.numberToLocation(i_segment); + i_segment = segment_set.segmentNumberToIndex(i_segment); m_segment_perforations_[i_segment].push_back(i); temp_perf_depth[i] = completion_set.get(i).getCenterDepth(); } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 23a872975..15de839d2 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -165,7 +166,7 @@ namespace Opm } } } - // perfPressures + // perfPressures if( num_perf_old_well == num_perf_this_well ) { int oldPerf_idx = oldPerf_idx_beg; @@ -206,6 +207,17 @@ namespace Opm } } } + + { + // we need to create a trival segment related values to avoid there will be some + // multi-segment wells added later. + top_segment_index_.reserve(nw); + for (int w = 0; w < nw; ++w) { + top_segment_index_.push_back(w); + } + segpress_ = bhp(); + segrates_ = wellRates(); + } } template @@ -247,7 +259,6 @@ namespace Opm if (pu.has_solvent) { // add solvent component for( int w = 0; w < nw; ++w ) { - using rt = data::Rates::opt; res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) ); } } @@ -282,6 +293,197 @@ namespace Opm } + /// init the MS well related. + template + void initWellStateMSWell(const Wells* wells, const std::vector& wells_ecl, + const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state) + { + // still using the order in wells + const int nw = wells->number_of_wells; + if (nw == 0) { + return; + } + + top_segment_index_.clear(); + top_segment_index_.reserve(nw); + segpress_.clear(); + segpress_.reserve(nw); + segrates_.clear(); + segrates_.reserve(nw * numPhases()); + + nseg_ = 0; + // in the init function, the well rates and perforation rates have been initialized or copied from prevState + // what we do here, is to set the segment rates and perforation rates + for (int w = 0; w < nw; ++w) { + const int nw_wells_ecl = wells_ecl.size(); + int index_well_ecl = 0; + const std::string well_name(wells->name[w]); + for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { + if (well_name == wells_ecl[index_well_ecl]->name()) { + break; + } + } + + // It should be able to find in wells_ecl. + if (index_well_ecl == 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_ecl]; + top_segment_index_.push_back(nseg_); + if ( !well_ecl->isMultiSegment(time_step) ) { // not multi-segment well + nseg_ += 1; + segpress_.push_back(bhp()[w]); + const int np = numPhases(); + for (int p = 0; p < np; ++p) { + segrates_.push_back(wellRates()[np * w + p]); + } + } else { // it is a multi-segment well + const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step); + // assuming the order of the perforations in well_ecl is the same with Wells + const CompletionSet& completion_set = well_ecl->getCompletions(time_step); + // number of segment for this single well + const int well_nseg = segment_set.numberSegment(); + const int nperf = completion_set.size(); + nseg_ += well_nseg; + // we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment + // that is why I think we should use a well model to initialize the WellState here + std::vector> segment_perforations(well_nseg); + for (int perf = 0; perf < nperf; ++perf) { + const Completion& completion = completion_set.get(perf); + const int segment_number = completion.getSegmentNumber(); + const int segment_index = segment_set.segmentNumberToIndex(segment_number); + segment_perforations[segment_index].push_back(perf); + } + + std::vector> segment_inlets(well_nseg); + for (int seg = 0; seg < well_nseg; ++seg) { + const Segment& segment = segment_set[seg]; + const int segment_number = segment.segmentNumber(); + const int outlet_segment_number = segment.outletSegment(); + if (outlet_segment_number > 0) { + const int segment_index = segment_set.segmentNumberToIndex(segment_number); + const int outlet_segment_index = segment_set.segmentNumberToIndex(outlet_segment_number); + segment_inlets[outlet_segment_index].push_back(segment_index); + } + } + + + // for the segrates_, now it becomes a recursive solution procedure. + { + const int np = numPhases(); + const int start_perf = wells->well_connpos[w]; + const int start_perf_next_well = wells->well_connpos[w + 1]; + assert(nperf == (start_perf_next_well - start_perf)); // make sure the information from wells_ecl consistent with wells + if (pu.phase_used[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + // scale the phase rates for Gas to avoid too bad initial guess for gas fraction + // it will probably benefit the standard well too, while it needs to be justified + // TODO: to see if this strategy can benefit StandardWell too + // TODO: it might cause big problem for gas rate control or if there is a gas rate limit + // maybe the best way is to initialize the fractions first then get the rates + for (int perf = 0; perf < nperf; perf++) { + const int perf_pos = start_perf + perf; + perfPhaseRates()[np * perf_pos + gaspos] *= 100.; + } + } + + const std::vector perforation_rates(perfPhaseRates().begin() + np * start_perf, + perfPhaseRates().begin() + np * start_perf_next_well); // the perforation rates for this well + std::vector segment_rates; + calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, segment_rates); + std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_)); + } + + // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment + // if there is no perforation associated with this segment, it uses the pressure from the outlet segment + // which requres the ordering is successful + // Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be + // improved during the solveWellEq process + { + // top segment is always the first one, and its pressure is the well bhp + segpress_.push_back(bhp()[w]); + const int top_segment = top_segment_index_[w]; + const int start_perf = wells->well_connpos[w]; + for (int seg = 1; seg < well_nseg; ++seg) { + if ( !segment_perforations[seg].empty() ) { + const int first_perf = segment_perforations[seg][0]; + segpress_.push_back(perfPress()[start_perf + first_perf]); + } else { + // segpress_.push_back(bhp); // may not be a good decision + // using the outlet segment pressure // it needs the ordering is correct + const int outlet_seg = segment_set[seg].outletSegment(); + segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]); + } + } + } + } + } + assert(int(segpress_.size()) == nseg_); + assert(int(segrates_.size()) == nseg_ * numPhases() ); + + if (!prev_well_state.wellMap().empty()) { + // copying MS well related + const auto& end = prev_well_state.wellMap().end(); + const int np = numPhases(); + for (int w = 0; w < nw; ++w) { + const std::string name( wells->name[w] ); + const auto& it = prev_well_state.wellMap().find( name ); + + if (it != end) { // the well is found in the prev_well_state + // TODO: the well with same name can change a lot, like they might not have same number of segments + // we need to handle that later. + // for now, we just copy them. + const int old_index_well = (*it).second[0]; + const int new_index_well = w; + const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well); + const int new_top_segmnet_index = topSegmentIndex(new_index_well); + int number_of_segment = 0; + // if it is the last well in list + if (new_index_well == int(top_segment_index_.size()) - 1) { + number_of_segment = nseg_ - new_top_segmnet_index; + } else { + number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index; + } + + for (int i = 0; i < number_of_segment * np; ++i) { + segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i]; + } + + for (int i = 0; i < number_of_segment; ++i) { + segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i]; + } + } + } + } + } + + + static void calculateSegmentRates(const std::vector>& segment_inlets, const std::vector>&segment_perforations, + const std::vector& perforation_rates, const int np, const int segment, std::vector& segment_rates) + { + // the rate of the segment equals to the sum of the contribution from the perforations and inlet segment rates. + // the first segment is always the top segment, its rates should be equal to the well rates. + assert(segment_inlets.size() == segment_perforations.size()); + const int well_nseg = segment_inlets.size(); + if (segment == 0) { // beginning the calculation + segment_rates.resize(np * well_nseg, 0.0); + } + // contributions from the perforations belong to this segment + for (const int& perf : segment_perforations[segment]) { + for (int p = 0; p < np; ++p) { + segment_rates[np * segment + p] += perforation_rates[np * perf + p]; + } + } + for (const int& inlet_seg : segment_inlets[segment]) { + calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates); + for (int p = 0; p < np; ++p) { + segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p]; + } + } + } + + bool isNewWell(const int w) const { return is_new_well_[w]; } @@ -305,6 +507,38 @@ namespace Opm return solvent_well_rate; } + const std::vector& segRates() const + { + return segrates_; + } + + std::vector& segRates() + { + return segrates_; + } + + const std::vector& segPress() const + { + return segpress_; + } + + std::vector& segPress() + { + return segpress_; + } + + int numSegment() const + { + return nseg_; + } + + int topSegmentIndex(const int w) const + { + assert(w < int(top_segment_index_.size()) ); + + return top_segment_index_[w]; + } + private: std::vector perfphaserates_; std::vector current_controls_; @@ -315,6 +549,16 @@ namespace Opm // will have very wrong compositions for production wells, will mostly cause // problem with VFP interpolation std::vector is_new_well_; + + // MS well related + // for StandardWell, the number of segments will be one + std::vector segrates_; + std::vector segpress_; + // the index of the top segments, which is used to locate the + // multisegment well related information in WellState + std::vector top_segment_index_; + int nseg_; // total number of the segments + }; } // namespace Opm diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index 77fab123a..a66e76270 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -47,6 +47,7 @@ #include #include +#include #include #include @@ -56,7 +57,7 @@ #include #include -#include +#include // maybe should just include BlackoilModelEbos.hpp namespace Ewoms { @@ -122,9 +123,10 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) { 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); + const Opm::BlackoilModelParameters param; + BOOST_CHECK_THROW( StandardWell( well, -1, wells, param), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param), std::invalid_argument); } @@ -137,6 +139,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { { const int nw = wells_struct ? (wells_struct->number_of_wells) : 0; + const Opm::BlackoilModelParameters param; for (int w = 0; w < nw; ++w) { const std::string well_name(wells_struct->name[w]); @@ -150,7 +153,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { // 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) ); + wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param) ); } }