diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 74d7a150c..76cf3f793 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -63,6 +63,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/GasLiftStage2.cpp opm/simulators/wells/GlobalWellInfo.cpp opm/simulators/wells/GroupState.cpp + opm/simulators/wells/MultisegmentWellEval.cpp opm/simulators/wells/MultisegmentWellGeneric.cpp opm/simulators/wells/ParallelWellInfo.cpp opm/simulators/wells/SegmentState.cpp diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index f41060e05..8258f998c 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -29,7 +29,10 @@ #include #include #include +#include +#include #include + #if HAVE_UMFPACK #include #endif // HAVE_UMFPACK diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 97958a5f2..59a7fe095 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -23,7 +23,7 @@ #define OPM_MULTISEGMENTWELL_HEADER_INCLUDED #include -#include +#include #include @@ -33,10 +33,15 @@ namespace Opm template class MultisegmentWell : public WellInterface - , public MultisegmentWellGeneric> + , public MultisegmentWellEval, + GetPropType, + GetPropType> { public: - typedef WellInterface Base; + using Base = WellInterface; + using MSWEval = MultisegmentWellEval, + GetPropType, + GetPropType>; using typename Base::Simulator; using typename Base::IntensiveQuantities; @@ -62,63 +67,21 @@ namespace Opm using Base::Oil; using Base::Gas; - // 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. - - // Table showing the primary variable indices, depending on what phases are present: - // - // WOG OG WG WO W/O/G (single phase) - // GTotal 0 0 0 0 0 - // WFrac 1 -1000 1 1 -1000 - // GFrac 2 1 -1000 -1000 -1000 - // Spres 3 2 2 2 1 - - static constexpr bool has_water = (Indices::waterSaturationIdx >= 0); - static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0); - static constexpr bool has_oil = (numPhases - has_gas - has_water) > 0; - - // In the implementation, one should use has_wfrac_variable - // rather than has_water to check if you should do something - // with the variable at the WFrac location, similar for GFrac. - static constexpr bool has_wfrac_variable = has_water && numPhases > 1; - static constexpr bool has_gfrac_variable = has_gas && has_oil; - - static constexpr int GTotal = 0; - static constexpr int WFrac = has_wfrac_variable ? 1 : -1000; - static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000; - static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1; - - // the number of well equations TODO: it should have a more general strategy for it - static const int numWellEq = numPhases + 1; - using typename Base::Scalar; /// the matrix and vector types for the reservoir 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; + using typename MSWEval::EvalWell; + using typename MSWEval::BVectorWell; + using typename MSWEval::DiagMatWell; + using typename MSWEval::OffDiagMatrixBlockWellType; + using MSWEval::GFrac; + using MSWEval::WFrac; + using MSWEval::GTotal; + using MSWEval::SPres; + using MSWEval::numWellEq; MultisegmentWell(const Well& well, const ParallelWellInfo& pw_info, @@ -137,7 +100,6 @@ namespace Opm const int num_cells, const std::vector< Scalar >& B_avg) override; - virtual void initPrimaryVariablesEvaluation() const override; virtual void gasLiftOptimizationStage1 ( @@ -157,18 +119,16 @@ namespace Opm DeferredLogger& deferred_logger) const; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, DeferredLogger& deferred_logger, const bool relax_tolerance = false) const override; + virtual ConvergenceReport getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + DeferredLogger& deferred_logger, + const bool relax_tolerance = false) const override; /// Ax = Ax - C D^-1 B x virtual void apply(const BVector& x, BVector& Ax) const override; /// r = r - C D^-1 Rw virtual void apply(BVector& r) const override; -#if HAVE_CUDA || HAVE_OPENCL - /// add the contribution (C, D, B matrices) of this Well to the WellContributions object - void addWellContribution(WellContributions& wellContribs) const; -#endif - /// using the solution x to recover the solution xw for wells and applying /// xw to update Well State virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, @@ -260,67 +220,11 @@ namespace Opm using Base::calculateBhpFromThp; using Base::getALQ; - // 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. It has offdiagonal entries for inlets and outlets. - mutable DiagMatWell duneD_; - /// \brief solver for diagonal matrix - /// - /// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials... - mutable std::shared_ptr > duneDSolver_; - - // 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_; - // the intial amount of fluids in each segment under surface condition std::vector > segment_fluid_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_; - - // the upwinding segment for each segment based on the flow direction - std::vector upwinding_segments_; - mutable int debug_cost_counter_ = 0; - std::vector> segment_phase_fractions_; - - std::vector> segment_phase_viscosities_; - - std::vector> segment_phase_densities_; - - - void initMatrixAndVectors(const int num_cells) const; - - EvalWell getBhp() const; - EvalWell getQs(const int comp_idx) const; - EvalWell getWQTotal() const; - - // xw = inv(D)*(rw - C*x) - void recoverSolutionWell(const BVector& x, BVectorWell& xw) const; - // updating the well_state based on well solution dwells void updateWellState(const BVectorWell& dwells, WellState& well_state, @@ -334,16 +238,6 @@ namespace Opm // 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 unsigned 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 computePerfRatePressure(const IntensiveQuantities& int_quants, const std::vector& mob_perfcells, const double Tw, @@ -357,23 +251,10 @@ namespace Opm double& perf_vap_oil_rate, DeferredLogger& deferred_logger) const; - // convert a Eval from reservoir to contain the derivative related to wells - EvalWell extendEval(const Eval& in) const; - - void updateThp(WellState& well_state, DeferredLogger& deferred_logger) 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 getSegmentRateUpwinding(const int seg, const size_t comp_idx) const; - - EvalWell getSegmentGTotal(const int seg) const; - // get the mobility for specific perforation void getMobility(const Simulator& ebosSimulator, const int perf, @@ -392,32 +273,6 @@ namespace Opm computeWellPotentialWithTHP(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const; - void assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::InjectionControls& inj_controls, - const Well::ProductionControls& prod_controls, - DeferredLogger& deferred_logger); - - void assemblePressureEq(const int seg, const UnitSystem& unit_system, - WellState& well_state, DeferredLogger& deferred_logger) const; - - void assembleDefaultPressureEq(const int seg, WellState& well_state) const; - - // hytrostatic pressure loss - EvalWell getHydroPressureLoss(const int seg) const; - - // frictinal pressure loss - EvalWell getFrictionPressureLoss(const int seg) const; - - void handleAccelerationPressureLoss(const int seg, WellState& well_state) const; - - // handling the overshooting and undershooting of the fractions - void processFractions(const int seg) const; - - void updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const; - virtual double getRefDensity() const override; virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator, @@ -440,21 +295,6 @@ namespace Opm EvalWell getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const; - std::vector getWellResiduals(const std::vector& B_avg, - DeferredLogger& deferred_logger) const; - - double getResidualMeasureValue(const WellState& well_state, - const std::vector& residuals, - DeferredLogger& deferred_logger) const; - - double getControlTolerance(const WellState& well_state, DeferredLogger& deferred_logger) const; - - void checkConvergenceControlEq(const WellState& well_state, - ConvergenceReport& report, - DeferredLogger& deferred_logger) const; - - void updateUpwindingSegments(); - // turn on crossflow to avoid singular well equations // when the well is banned from cross-flow and the BHP is not properly initialized, // we turn on crossflow to avoid singular well equations. It can result in wrong-signed @@ -477,19 +317,6 @@ namespace Opm double maxPerfPress(const Simulator& ebos_simulator) const; - // pressure drop for Spiral ICD segment (WSEGSICD) - EvalWell pressureDropSpiralICD(const int seg) const; - - // pressure drop for Autonomous ICD segment (WSEGAICD) - EvalWell pressureDropAutoICD(const int seg, const UnitSystem& unit_system) const; - - // pressure drop for sub-critical valve (WSEGVALV) - EvalWell pressureDropValve(const int seg) const; - - // assemble pressure equation for ICD segments - void assembleICDPressureEq(const int seg, const UnitSystem& unit_system, - WellState& well_state, DeferredLogger& deferred_logger) const; - // check whether the well is operable under BHP limit with current reservoir condition virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) override; @@ -498,7 +325,6 @@ namespace Opm // updating the inflow based on the current reservoir condition virtual void updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const override; - }; } diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp new file mode 100644 index 000000000..3155508e9 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -0,0 +1,2017 @@ +/* + 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 +#include + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace Opm +{ + +template +MultisegmentWellEval:: +MultisegmentWellEval(WellInterfaceIndices& baseif) + : MultisegmentWellGeneric(baseif) + , baseif_(baseif) + , upwinding_segments_(this->numberOfSegments(), 0) + , segment_densities_(this->numberOfSegments(), 0.0) + , segment_mass_rates_(this->numberOfSegments(), 0.0) + , segment_viscosities_(this->numberOfSegments(), 0.0) + , segment_phase_densities_(this->numberOfSegments(), std::vector(baseif_.numComponents(), 0.0)) // number of phase here? + , segment_phase_fractions_(this->numberOfSegments(), std::vector(baseif_.numComponents(), 0.0)) // number of phase here? + , segment_phase_viscosities_(this->numberOfSegments(), std::vector(baseif_.numComponents(), 0.0)) // number of phase here? + , cell_perforation_depth_diffs_(baseif_.numPerfs(), 0.0) + , cell_perforation_pressure_diffs_(baseif_.numPerfs(), 0.0) +{ +} + +template +void +MultisegmentWellEval:: +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 = this->numberOfSegments(); + for (const std::vector& inlets : this->segment_inlets_) { + nnz_d += 2 * inlets.size(); + } + duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d); + } + duneB_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); + duneC_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); + + // 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 = this->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 = this->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 : this->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 : this->segment_perforations_[row.index()]) { + const int cell_idx = baseif_.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 : this->segment_perforations_[row.index()]) { + const int cell_idx = baseif_.cells()[perf]; + row.insert(cell_idx); + } + } + + resWell_.resize(this->numberOfSegments()); + + primary_variables_.resize(this->numberOfSegments()); + primary_variables_evaluation_.resize(this->numberOfSegments()); +} + +template +void +MultisegmentWellEval:: +initPrimaryVariablesEvaluation() const +{ + for (int seg = 0; seg < this->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 + Indices::numEq, 1.0); + } + } +} + +template +void +MultisegmentWellEval:: +checkConvergenceControlEq(const WellState& well_state, + ConvergenceReport& report, + const double tolerance_pressure_ms_wells, + const double tolerance_wells, + const double max_residual_allowed, + DeferredLogger& deferred_logger) const +{ + double control_tolerance = 0.; + using CR = ConvergenceReport; + CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid; + + const int well_index = baseif_.indexOfWell(); + if (baseif_.isInjector() ) + { + auto current = well_state.currentInjectionControl(well_index); + switch(current) { + case Well::InjectorCMode::THP: + ctrltype = CR::WellFailure::Type::ControlTHP; + control_tolerance = tolerance_pressure_ms_wells; + break; + case Well::InjectorCMode::BHP: + ctrltype = CR::WellFailure::Type::ControlBHP; + control_tolerance = tolerance_pressure_ms_wells; + break; + case Well::InjectorCMode::RATE: + case Well::InjectorCMode::RESV: + ctrltype = CR::WellFailure::Type::ControlRate; + control_tolerance = tolerance_wells; + break; + case Well::InjectorCMode::GRUP: + ctrltype = CR::WellFailure::Type::ControlRate; + control_tolerance = tolerance_wells; + break; + default: + OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger); + } + } + + if (baseif_.isProducer() ) + { + auto current = well_state.currentProductionControl(well_index); + switch(current) { + case Well::ProducerCMode::THP: + ctrltype = CR::WellFailure::Type::ControlTHP; + control_tolerance = tolerance_pressure_ms_wells; + break; + case Well::ProducerCMode::BHP: + ctrltype = CR::WellFailure::Type::ControlBHP; + control_tolerance = tolerance_pressure_ms_wells; + break; + case Well::ProducerCMode::ORAT: + case Well::ProducerCMode::WRAT: + case Well::ProducerCMode::GRAT: + case Well::ProducerCMode::LRAT: + case Well::ProducerCMode::RESV: + case Well::ProducerCMode::CRAT: + ctrltype = CR::WellFailure::Type::ControlRate; + control_tolerance = tolerance_wells; + break; + case Well::ProducerCMode::GRUP: + ctrltype = CR::WellFailure::Type::ControlRate; + control_tolerance = tolerance_wells; + break; + default: + OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger); + } + } + + const double well_control_residual = std::abs(resWell_[0][SPres]); + const int dummy_component = -1; + if (std::isnan(well_control_residual)) { + report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()}); + } else if (well_control_residual > max_residual_allowed * 10.) { + report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, baseif_.name()}); + } else if ( well_control_residual > control_tolerance) { + report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, baseif_.name()}); + } +} + +template +ConvergenceReport +MultisegmentWellEval:: +getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + DeferredLogger& deferred_logger, + const double max_residual_allowed, + const double tolerance_wells, + const double relaxed_inner_tolerance_flow_ms_well, + const double tolerance_pressure_ms_wells, + const double relaxed_inner_tolerance_pressure_ms_well, + const bool relax_tolerance) const +{ + assert(int(B_avg.size()) == baseif_.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(this->numberOfSegments(), std::vector(numWellEq, 0.0)); + for (int seg = 0; seg < this->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 eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + if (eq_idx < baseif_.numComponents()) { // phase or component mass equations + const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx]; + if (flux_residual > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = flux_residual; + } + } else { // pressure or control equation + // for the top segment (seg == 0), it is control equation, will be checked later separately + if (seg > 0) { + // Pressure equation + const double pressure_residual = abs_residual[seg][eq_idx]; + if (pressure_residual > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = pressure_residual; + } + } + } + } + } + + using CR = ConvergenceReport; + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + if (eq_idx < baseif_.numComponents()) { // phase or component mass equations + const double flux_residual = maximum_residual[eq_idx]; + // TODO: the report can not handle the segment number yet. + + if (std::isnan(flux_residual)) { + report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, baseif_.name()}); + } else if (flux_residual > max_residual_allowed) { + report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, baseif_.name()}); + } else if (!relax_tolerance && flux_residual > tolerance_wells) { + report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, baseif_.name()}); + } else if (flux_residual > relaxed_inner_tolerance_flow_ms_well) { + report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, baseif_.name()}); + } + } else { // pressure equation + const double pressure_residual = maximum_residual[eq_idx]; + const int dummy_component = -1; + if (std::isnan(pressure_residual)) { + report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, baseif_.name()}); + } else if (std::isinf(pressure_residual)) { + report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, baseif_.name()}); + } else if (!relax_tolerance && pressure_residual > tolerance_pressure_ms_wells) { + report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, baseif_.name()}); + } else if (pressure_residual > relaxed_inner_tolerance_pressure_ms_well) { + report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, baseif_.name()}); + } + } + } + + checkConvergenceControlEq(well_state, + report, + tolerance_pressure_ms_wells, + tolerance_wells, + max_residual_allowed, + deferred_logger); + + return report; +} + +template +void +MultisegmentWellEval:: +processFractions(const int seg) const +{ + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Gas = BlackoilPhases::Vapour; + + const PhaseUsage& pu = baseif_.phaseUsage(); + + std::vector fractions(baseif_.numPhases(), 0.0); + + assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); + const int oil_pos = pu.phase_pos[Oil]; + fractions[oil_pos] = 1.0; + + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int water_pos = pu.phase_pos[Water]; + fractions[water_pos] = primary_variables_[seg][WFrac]; + fractions[oil_pos] -= fractions[water_pos]; + } + + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int gas_pos = pu.phase_pos[Gas]; + fractions[gas_pos] = primary_variables_[seg][GFrac]; + fractions[oil_pos] -= fractions[gas_pos]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int water_pos = pu.phase_pos[Water]; + if (fractions[water_pos] < 0.0) { + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]); + } + fractions[oil_pos] /= (1.0 - fractions[water_pos]); + fractions[water_pos] = 0.0; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int gas_pos = pu.phase_pos[Gas]; + if (fractions[gas_pos] < 0.0) { + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + 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 ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]); + } + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]); + } + fractions[oil_pos] = 0.0; + } + + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]]; + } + + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]]; + } +} + +template +void +MultisegmentWellEval:: +updateWellState(const BVectorWell& dwells, + const double relaxation_factor, + const double dFLimit, + const double max_pressure_change) const +{ + const std::vector > old_primary_variables = primary_variables_; + + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + if (has_wfrac_variable) { + 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 (has_gfrac_variable) { + 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] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); + } + + // 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]; + + // make sure that no injector produce and no producer inject + if (seg == 0) { + if (baseif_.isInjector()) { + primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0); + } else { + primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0); + } + } + } + } +} + +template +void +MultisegmentWellEval:: +updatePrimaryVariables(const WellState& well_state) const +{ + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Gas = BlackoilPhases::Vapour; + + // TODO: to test using rate conversion coefficients to see if it will be better than + // this default one + if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return; + + const Well& well = baseif_.wellEcl(); + + // the index of the top segment in the WellState + const auto& segments = well_state.segments(baseif_.indexOfWell()); + const auto& segment_rates = segments.rates; + const auto& segment_pressure = segments.pressure; + const PhaseUsage& pu = baseif_.phaseUsage(); + + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + // calculate the total rate for each segment + double total_seg_rate = 0.0; + // the segment pressure + primary_variables_[seg][SPres] = segment_pressure[seg]; + // 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 < baseif_.numPhases(); p++) { + total_seg_rate += baseif_.scalingFactor(p) * segment_rates[baseif_.numPhases() * seg + p]; + } + + primary_variables_[seg][GTotal] = total_seg_rate; + if (std::abs(total_seg_rate) > 0.) { + if (has_wfrac_variable) { + const int water_pos = pu.phase_pos[Water]; + primary_variables_[seg][WFrac] = baseif_.scalingFactor(water_pos) * segment_rates[baseif_.numPhases() * seg + water_pos] / total_seg_rate; + } + if (has_gfrac_variable) { + const int gas_pos = pu.phase_pos[Gas]; + primary_variables_[seg][GFrac] = baseif_.scalingFactor(gas_pos) * segment_rates[baseif_.numPhases() * seg + gas_pos] / total_seg_rate; + } + } else { // total_seg_rate == 0 + if (baseif_.isInjector()) { + // only single phase injection handled + auto phase = well.getInjectionProperties().injectorType; + + if (has_wfrac_variable) { + if (phase == InjectorType::WATER) { + primary_variables_[seg][WFrac] = 1.0; + } else { + primary_variables_[seg][WFrac] = 0.0; + } + } + + if (has_gfrac_variable) { + if (phase == InjectorType::GAS) { + primary_variables_[seg][GFrac] = 1.0; + } else { + primary_variables_[seg][GFrac] = 0.0; + } + } + + } else if (baseif_.isProducer()) { // producers + if (has_wfrac_variable) { + primary_variables_[seg][WFrac] = 1.0 / baseif_.numPhases(); + } + + if (has_gfrac_variable) { + primary_variables_[seg][GFrac] = 1.0 / baseif_.numPhases(); + } + } + } + } +} + +template +void +MultisegmentWellEval:: +recoverSolutionWell(const BVector& x, BVectorWell& xw) const +{ + if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return; + + BVectorWell resWell = resWell_; + // resWell = resWell - B * x + duneB_.mmv(x, resWell); + // xw = D^-1 * resWell + xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +volumeFraction(const int seg, + const unsigned compIdx) const +{ + if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { + return primary_variables_evaluation_[seg][WFrac]; + } + + if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { + return primary_variables_evaluation_[seg][GFrac]; + } + + // Oil fraction + EvalWell oil_fraction = 1.0; + if (has_wfrac_variable) { + oil_fraction -= primary_variables_evaluation_[seg][WFrac]; + } + + if (has_gfrac_variable) { + oil_fraction -= primary_variables_evaluation_[seg][GFrac]; + } + /* if (has_solvent) { + oil_fraction -= primary_variables_evaluation_[seg][SFrac]; + } */ + return oil_fraction; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +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 = baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + if (scale > 0.) { + return volumeFraction(seg, comp_idx) / scale; + } + + return volumeFraction(seg, comp_idx); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +surfaceVolumeFraction(const int seg, + const int comp_idx) const +{ + EvalWell sum_volume_fraction_scaled = 0.; + for (int idx = 0; idx < baseif_.numComponents(); ++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 +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getSegmentRateUpwinding(const int seg, + const size_t comp_idx) const +{ + const int seg_upwind = upwinding_segments_[seg]; + // the result will contain the derivative with resepct to GTotal in segment seg, + // and the derivatives with respect to WFrac GFrac in segment seg_upwind. + // the derivative with respect to SPres should be zero. + if (seg == 0 && baseif_.isInjector()) { + const Well& well = baseif_.wellEcl(); + auto phase = well.getInjectionProperties().injectorType; + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) + && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx + && phase == InjectorType::WATER) + return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) + && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx + && phase == InjectorType::OIL) + return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) + && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx + && phase == InjectorType::GAS) + return primary_variables_evaluation_[seg][GTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + + return 0.0; + } + + const EvalWell segment_rate = primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx); + + assert(segment_rate.derivative(SPres + Indices::numEq) == 0.); + + return segment_rate; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +extendEval(const Eval& in) const +{ + EvalWell out = 0.0; + out.setValue(in.value()); + for(int eq_idx = 0; eq_idx < Indices::numEq;++eq_idx) { + out.setDerivative(eq_idx, in.derivative(eq_idx)); + } + return out; +} + +template +void +MultisegmentWellEval:: +computeSegmentFluidProperties(const EvalWell& temperature, + const EvalWell& saltConcentration, + int pvt_region_index) +{ + std::vector surf_dens(baseif_.numComponents()); + // Surface density. + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index); + } + + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + // the compostion of the components inside wellbore under surface condition + std::vector mix_s(baseif_.numComponents(), 0.0); + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); + } + + std::vector b(baseif_.numComponents(), 0.0); + std::vector visc(baseif_.numComponents(), 0.0); + std::vector& phase_densities = segment_phase_densities_[seg]; + + const EvalWell seg_pressure = getSegmentPressure(seg); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b[waterCompIdx] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration); + visc[waterCompIdx] = + FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration); + // TODO: double check here + // TODO: should not we use phaseIndex here? + phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx]; + } + + EvalWell rv(0.0); + // gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[oilCompIdx] > 0.0) { + if (mix_s[gasCompIdx] > 0.0) { + rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; + } + + if (rv > rvmax) { + rv = rvmax; + } + b[gasCompIdx] = + FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); + visc[gasCompIdx] = + FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv); + phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx] + + rv * b[gasCompIdx] * surf_dens[oilCompIdx]; + } else { // no oil exists + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[gasCompIdx] = + FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx]; + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[gasCompIdx] = + FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } + + EvalWell rs(0.0); + // oil phase + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[gasCompIdx] > 0.0) { + if (mix_s[oilCompIdx] > 0.0) { + rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; + } + + if (rs > rsmax) { + rs = rsmax; + } + b[oilCompIdx] = + FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); + visc[oilCompIdx] = + FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs); + phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx] + + rs * b[oilCompIdx] * surf_dens[gasCompIdx]; + } else { // no oil exists + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[oilCompIdx] = + FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx]; + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + visc[oilCompIdx] = + FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); + } + } + + segment_phase_viscosities_[seg] = visc; + + std::vector mix(mix_s); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + + const EvalWell d = 1.0 - rs * rv; + + if (rs != 0.0) { // rs > 0.0? + mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; + } + if (rv != 0.0) { // rv > 0.0? + mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; + } + } + + EvalWell volrat(0.0); + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + volrat += mix[comp_idx] / b[comp_idx]; + } + + this->segment_viscosities_[seg] = 0.; + // calculate the average viscosity + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; + // TODO: a little more work needs to be done to handle the negative fractions here + this->segment_phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; + this->segment_viscosities_[seg] += visc[comp_idx] * this->segment_phase_fractions_[seg][comp_idx]; + } + + EvalWell density(0.0); + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + density += surf_dens[comp_idx] * mix_s[comp_idx]; + } + this->segment_densities_[seg] = density / volrat; + + // calculate the mass rates + segment_mass_rates_[seg] = 0.; + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + const EvalWell rate = getSegmentRateUpwinding(seg, comp_idx); + this->segment_mass_rates_[seg] += rate * surf_dens[comp_idx]; + } + } +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getSegmentPressure(const int seg) const +{ + return primary_variables_evaluation_[seg][SPres]; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getBhp() const +{ + return getSegmentPressure(0); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getSegmentRate(const int seg, + const int comp_idx) const +{ + return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getQs(const int comp_idx) const +{ + return getSegmentRate(0, comp_idx); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getSegmentGTotal(const int seg) const +{ + return primary_variables_evaluation_[seg][GTotal]; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getWQTotal() const +{ + return getSegmentGTotal(0); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getHydroPressureLoss(const int seg) const +{ + return segment_densities_[seg] * baseif_.gravity() * this->segment_depth_diffs_[seg]; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getFrictionPressureLoss(const int seg) const +{ + const EvalWell mass_rate = segment_mass_rates_[seg]; + const int seg_upwind = upwinding_segments_[seg]; + EvalWell density = segment_densities_[seg_upwind]; + EvalWell visc = segment_viscosities_[seg_upwind]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (seg != seg_upwind) { + density.clearDerivatives(); + visc.clearDerivatives(); + } + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); + const double length = this->segmentSet()[seg].totalLength() - this->segmentSet()[outlet_segment_index].totalLength(); + assert(length > 0.); + const double roughness = this->segmentSet()[seg].roughness(); + const double area = this->segmentSet()[seg].crossArea(); + const double diameter = this->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 +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +pressureDropSpiralICD(const int seg) const +{ + const SICD& sicd = this->segmentSet()[seg].spiralICD(); + + const int seg_upwind = upwinding_segments_[seg]; + const std::vector& phase_fractions = segment_phase_fractions_[seg_upwind]; + const std::vector& phase_viscosities = segment_phase_viscosities_[seg_upwind]; + + EvalWell water_fraction = 0.; + EvalWell water_viscosity = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + water_fraction = phase_fractions[water_pos]; + water_viscosity = phase_viscosities[water_pos]; + } + + EvalWell oil_fraction = 0.; + EvalWell oil_viscosity = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + oil_fraction = phase_fractions[oil_pos]; + oil_viscosity = phase_viscosities[oil_pos]; + } + + EvalWell gas_fraction = 0.; + EvalWell gas_viscosity = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + gas_fraction = phase_fractions[gas_pos]; + gas_viscosity = phase_viscosities[gas_pos]; + } + + EvalWell density = segment_densities_[seg_upwind]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (seg != seg_upwind) { + water_fraction.clearDerivatives(); + water_viscosity.clearDerivatives(); + oil_fraction.clearDerivatives(); + oil_viscosity.clearDerivatives(); + gas_fraction.clearDerivatives(); + gas_viscosity.clearDerivatives(); + density.clearDerivatives(); + } + + const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, + oil_fraction, oil_viscosity, sicd); + const EvalWell mixture_viscosity = (water_fraction + oil_fraction) * liquid_emulsion_viscosity + gas_fraction * gas_viscosity; + + const EvalWell reservoir_rate = segment_mass_rates_[seg] / density; + + const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor(); + + const double viscosity_cali = sicd.viscosityCalibration(); + + using MathTool = MathToolbox; + + const double density_cali = sicd.densityCalibration(); + const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75); + const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25); + + // formulation before 2016, base_strength is used + // const double base_strength = sicd.strength() / density_cali; + // formulation since 2016, strength is used instead + const double strength = sicd.strength(); + + const double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0; + + return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +pressureDropAutoICD(const int seg, + const UnitSystem& unit_system) const +{ + const AutoICD& aicd = this->segmentSet()[seg].autoICD(); + + const int seg_upwind = this->upwinding_segments_[seg]; + const std::vector& phase_fractions = this->segment_phase_fractions_[seg_upwind]; + const std::vector& phase_viscosities = this->segment_phase_viscosities_[seg_upwind]; + const std::vector& phase_densities = this->segment_phase_densities_[seg_upwind]; + + EvalWell water_fraction = 0.; + EvalWell water_viscosity = 0.; + EvalWell water_density = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + water_fraction = phase_fractions[water_pos]; + water_viscosity = phase_viscosities[water_pos]; + water_density = phase_densities[water_pos]; + } + + EvalWell oil_fraction = 0.; + EvalWell oil_viscosity = 0.; + EvalWell oil_density = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + oil_fraction = phase_fractions[oil_pos]; + oil_viscosity = phase_viscosities[oil_pos]; + oil_density = phase_densities[oil_pos]; + } + + EvalWell gas_fraction = 0.; + EvalWell gas_viscosity = 0.; + EvalWell gas_density = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + gas_fraction = phase_fractions[gas_pos]; + gas_viscosity = phase_viscosities[gas_pos]; + gas_density = phase_densities[gas_pos]; + } + + EvalWell density = segment_densities_[seg_upwind]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (seg != seg_upwind) { + water_fraction.clearDerivatives(); + water_viscosity.clearDerivatives(); + water_density.clearDerivatives(); + oil_fraction.clearDerivatives(); + oil_viscosity.clearDerivatives(); + oil_density.clearDerivatives(); + gas_fraction.clearDerivatives(); + gas_viscosity.clearDerivatives(); + gas_density.clearDerivatives(); + density.clearDerivatives(); + } + + using MathTool = MathToolbox; + const EvalWell mixture_viscosity = MathTool::pow(water_fraction, aicd.waterViscExponent()) * water_viscosity + + MathTool::pow(oil_fraction, aicd.oilViscExponent()) * oil_viscosity + + MathTool::pow(gas_fraction, aicd.gasViscExponent()) * gas_viscosity; + + const EvalWell mixture_density = MathTool::pow(water_fraction, aicd.waterDensityExponent()) * water_density + + MathTool::pow(oil_fraction, aicd.oilDensityExponent()) * oil_density + + MathTool::pow(gas_fraction, aicd.gasDensityExponent()) * gas_density; + + const double rho_reference = aicd.densityCalibration(); + const double visc_reference = aicd.viscosityCalibration(); + const auto volume_rate_icd = this->segment_mass_rates_[seg] * aicd.scalingFactor() / mixture_density; + const double sign = volume_rate_icd <= 0. ? 1.0 : -1.0; + // convert 1 unit volume rate + using M = UnitSystem::measure; + const double unit_volume_rate = unit_system.to_si(M::geometric_volume_rate, 1.); + + // TODO: we did not consider the maximum allowed rate here + const auto result = sign / rho_reference * mixture_density * mixture_density + * MathTool::pow(visc_reference/mixture_viscosity, aicd.viscExponent()) + * aicd.strength() * MathTool::pow( -sign * volume_rate_icd, aicd.flowRateExponent()) + * std::pow(unit_volume_rate, (2. - aicd.flowRateExponent())) ; + return result; +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +pressureDropValve(const int seg) const +{ + const Valve& valve = this->segmentSet()[seg].valve(); + + const EvalWell& mass_rate = segment_mass_rates_[seg]; + const int seg_upwind = upwinding_segments_[seg]; + EvalWell visc = segment_viscosities_[seg_upwind]; + EvalWell density = segment_densities_[seg_upwind]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (seg != seg_upwind) { + visc.clearDerivatives(); + density.clearDerivatives(); + } + + const double additional_length = valve.pipeAdditionalLength(); + const double roughness = valve.pipeRoughness(); + const double diameter = valve.pipeDiameter(); + const double area = valve.pipeCrossArea(); + + const EvalWell friction_pressure_loss = + mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc); + + const double area_con = valve.conCrossArea(); + const double cv = valve.conFlowCoefficient(); + + const EvalWell constriction_pressure_loss = + mswellhelpers::valveContrictionPressureLoss(mass_rate, density, area_con, cv); + + const double sign = mass_rate <= 0. ? 1.0 : -1.0; + return sign * (friction_pressure_loss + constriction_pressure_loss); +} + +template +typename MultisegmentWellEval::EvalWell +MultisegmentWellEval:: +getSegmentSurfaceVolume(const EvalWell& temperature, + const EvalWell& saltConcentration, + const int pvt_region_index, + const int seg_idx) const +{ + const EvalWell seg_pressure = getSegmentPressure(seg_idx); + + std::vector mix_s(baseif_.numComponents(), 0.0); + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + mix_s[comp_idx] = surfaceVolumeFraction(seg_idx, comp_idx); + } + + std::vector b(baseif_.numComponents(), 0.); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b[waterCompIdx] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure, + saltConcentration); + } + + EvalWell rv(0.0); + // gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, + temperature, + seg_pressure); + if (rvmax < 0.0) { // negative rvmax can happen if the seg_pressure is outside the range of the table + rvmax = 0.0; + } + if (mix_s[oilCompIdx] > 0.0) { + if (mix_s[gasCompIdx] > 0.0) { + rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; + } + + if (rv > rvmax) { + rv = rvmax; + } + b[gasCompIdx] = + FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure, + rv); + } else { // no oil exists + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure); + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure); + } + } + + EvalWell rs(0.0); + // oil phase + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, + temperature, + seg_pressure); + if (rsmax < 0.0) { // negative rsmax can happen if the seg_pressure is outside the range of the table + rsmax = 0.0; + } + if (mix_s[gasCompIdx] > 0.0) { + if (mix_s[oilCompIdx] > 0.0) { + rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; + } + // std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl; + + if (rs > rsmax) { + rs = rsmax; + } + b[oilCompIdx] = + FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure, + rs); + } else { // no oil exists + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure); + } + } else { // no gas phase + // it is the same with zero mix_s[Gas] + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, + temperature, + seg_pressure); + } + } + + std::vector mix(mix_s); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + + const EvalWell d = 1.0 - rs * rv; + if (d <= 0.0 || d > 1.0) { + std::ostringstream sstr; + sstr << "Problematic d value " << d << " obtained for well " << baseif_.name() + << " during conversion to surface volume with rs " << rs + << ", rv " << rv << " and pressure " << seg_pressure + << " obtaining d " << d; + OpmLog::debug(sstr.str()); + OPM_THROW_NOLOG(NumericalIssue, sstr.str()); + } + + if (rs > 0.0) { // rs > 0.0? + mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; + } + if (rv > 0.0) { // rv > 0.0? + mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; + } + } + + EvalWell vol_ratio(0.0); + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + vol_ratio += mix[comp_idx] / b[comp_idx]; + } + + // We increase the segment volume with a factor 10 to stabilize the system. + const double volume = this->segmentSet()[seg_idx].volume(); + + return volume / vol_ratio; +} + +template +void +MultisegmentWellEval:: +computePerfRatePressure(const EvalWell& pressure_cell, + const EvalWell& rs, + const EvalWell& rv, + const std::vector& b_perfcells, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const +{ + std::vector cmix_s(baseif_.numComponents(), 0.0); + + // the composition of the components inside wellbore + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); + } + + + // pressure difference between the segment and the perforation + const EvalWell perf_seg_press_diff = baseif_.gravity() * segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf]; + // pressure difference between the perforation and the grid cell + const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; + + perf_press = pressure_cell - cell_perf_press_diff; + // 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 = perf_press - (segment_pressure + perf_seg_press_diff); + + // producing perforations + if ( drawdown > 0.0) { + // Do nothing is crossflow is not allowed + if (!allow_cf && baseif_.isInjector()) { + return; + } + + // compute component volumetric rates at standard conditions + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + const EvalWell cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); + cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const EvalWell cq_s_oil = cq_s[oilCompIdx]; + const EvalWell cq_s_gas = cq_s[gasCompIdx]; + cq_s[gasCompIdx] += rs * cq_s_oil; + cq_s[oilCompIdx] += rv * cq_s_gas; + } + } else { // injecting perforations + // Do nothing if crossflow is not allowed + if (!allow_cf && baseif_.isProducer()) { + return; + } + + // for injecting perforations, we use total mobility + EvalWell total_mob = mob_perfcells[0]; + for (int comp_idx = 1; comp_idx < baseif_.numComponents(); ++comp_idx) { + total_mob += mob_perfcells[comp_idx]; + } + + // injection perforations total volume rates + const EvalWell cqt_i = - Tw * (total_mob * drawdown); + + // compute volume ratio between connection and at standard conditions + EvalWell volume_ratio = 0.0; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // 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_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << baseif_.name() + << " during flux calculation" + << " with rs " << rs << " and rv " << rv, deferred_logger); + } + + const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; + volume_ratio += tmp_oil / b_perfcells[oilCompIdx]; + + const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d; + volume_ratio += tmp_gas / b_perfcells[gasCompIdx]; + } else { // not having gas and oil at the same time + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx]; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx]; + } + } + // injecting connections total volumerates at standard conditions + EvalWell cqt_is = cqt_i / volume_ratio; + for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { + cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; + } + } // end for injection perforations + + // calculating the perforation solution gas rate and solution oil rates + if (baseif_.isProducer()) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells + // s means standard condition, r means reservoir condition + // q_os = q_or * b_o + rv * q_gr * b_g + // q_gs = q_gr * g_g + rs * q_or * b_o + // d = 1.0 - rs * rv + // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) + // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) + + const double d = 1.0 - rv.value() * rs.value(); + // vaporized oil into gas + // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d + perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d; + // dissolved of gas in oil + // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d + perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d; + } + } +} + +template +void +MultisegmentWellEval:: +assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + const double rho, + DeferredLogger& deferred_logger) +{ + static constexpr int Gas = BlackoilPhases::Vapour; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Water = BlackoilPhases::Aqua; + + EvalWell control_eq(0.0); + + const auto& well = baseif_.wellEcl(); + + auto getRates = [&]() { + std::vector rates(3, 0.0); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)); + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); + } + return rates; + }; + + if (baseif_.wellIsStopped()) { + control_eq = getWQTotal(); + } else if (baseif_.isInjector() ) { + // Find scaling factor to get injection rate, + const InjectorType injectorType = inj_controls.injector_type; + double scaling = 1.0; + const auto& pu = baseif_.phaseUsage(); + switch (injectorType) { + case InjectorType::WATER: + { + scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]); + break; + } + case InjectorType::OIL: + { + scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]); + break; + } + case InjectorType::GAS: + { + scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]); + break; + } + default: + throw("Expected WATER, OIL or GAS as type for injectors " + well.name()); + } + const EvalWell injection_rate = getWQTotal() / scaling; + // Setup function for evaluation of BHP from THP (used only if needed). + auto bhp_from_thp = [&]() { + const auto rates = getRates(); + return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, rho, deferred_logger); + }; + // Call generic implementation. + baseif_.assembleControlEqInj(well_state, + group_state, + schedule, + summaryState, + inj_controls, + getBhp(), + injection_rate, + bhp_from_thp, + control_eq, + deferred_logger); + } else { + // Find rates. + const auto rates = getRates(); + // Setup function for evaluation of BHP from THP (used only if needed). + auto bhp_from_thp = [&]() { + return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, rho, deferred_logger); + }; + // Call generic implementation. + baseif_.assembleControlEqProd(well_state, + group_state, + schedule, + summaryState, + prod_controls, + getBhp(), + rates, + bhp_from_thp, + control_eq, + deferred_logger); + } + + // 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 + Indices::numEq); + } +} + +template +void +MultisegmentWellEval:: +updateThp(WellState& well_state, + const double rho, + DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = BlackoilPhases::Vapour; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Water = BlackoilPhases::Aqua; + + // When there is no vaild VFP table provided, we set the thp to be zero. + if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) { + well_state.update_thp(baseif_.indexOfWell(), 0.); + return; + } + + // the well is under other control types, we calculate the thp based on bhp and rates + std::vector rates(3, 0.0); + + const PhaseUsage& pu = baseif_.phaseUsage(); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + rates[ Water ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Water ] ]; + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + rates[ Oil ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Oil ] ]; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ]; + } + + const double bhp = well_state.bhp(baseif_.indexOfWell()); + + well_state.update_thp(baseif_.indexOfWell(), this->calculateThpFromBhp(rates, bhp, rho, deferred_logger)); +} + +template +void +MultisegmentWellEval:: +handleAccelerationPressureLoss(const int seg, + WellState& well_state) const +{ + const double area = this->segmentSet()[seg].crossArea(); + const EvalWell mass_rate = segment_mass_rates_[seg]; + const int seg_upwind = upwinding_segments_[seg]; + EvalWell density = segment_densities_[seg_upwind]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (seg != seg_upwind) { + density.clearDerivatives(); + } + + EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density); + // handling the velocity head of intlet segments + for (const int inlet : this->segment_inlets_[seg]) { + const int seg_upwind_inlet = upwinding_segments_[inlet]; + const double inlet_area = this->segmentSet()[inlet].crossArea(); + EvalWell inlet_density = this->segment_densities_[seg_upwind_inlet]; + // WARNING + // We disregard the derivatives from the upwind density to make sure derivatives + // wrt. to different segments dont get mixed. + if (inlet != seg_upwind_inlet) { + inlet_density.clearDerivatives(); + } + const EvalWell inlet_mass_rate = segment_mass_rates_[inlet]; + accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density); + } + + // We change the sign of the accelerationPressureLoss for injectors. + // Is this correct? Testing indicates that this is what the reference simulator does + const double sign = mass_rate < 0. ? 1.0 : - 1.0; + accelerationPressureLoss *= sign; + + well_state.segments(baseif_.indexOfWell()).pressure_drop_accel[seg] = accelerationPressureLoss.value(); + + resWell_[seg][SPres] -= accelerationPressureLoss.value(); + duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); + duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + Indices::numEq); + if (has_wfrac_variable) { + duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); + } + if (has_gfrac_variable) { + duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); + } +} + +template +void +MultisegmentWellEval:: +assembleDefaultPressureEq(const int seg, + WellState& well_state) 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 + // TODO: we might be able to add member variables to store these values, then we update well state + // after converged + const auto hydro_pressure_drop = getHydroPressureLoss(seg); + auto& segments = well_state.segments(baseif_.indexOfWell()); + segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value(); + pressure_equation -= hydro_pressure_drop; + + if (this->frictionalPressureLossConsidered()) { + const auto friction_pressure_drop = getFrictionPressureLoss(seg); + pressure_equation -= friction_pressure_drop; + segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); + } + + resWell_[seg][SPres] = pressure_equation.value(); + const int seg_upwind = upwinding_segments_[seg]; + duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); + duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + Indices::numEq); + if (has_wfrac_variable) { + duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); + } + if (has_gfrac_variable) { + duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); + } + + // contribution from the outlet segment + const int outlet_segment_index = this->segmentNumberToIndex(this->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 + Indices::numEq); + } + + if (this->accelerationalPressureLossConsidered()) { + handleAccelerationPressureLoss(seg, well_state); + } +} + +template +void +MultisegmentWellEval:: +updateWellStateFromPrimaryVariables(WellState& well_state, + const double rho, + DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = BlackoilPhases::Vapour; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Water = BlackoilPhases::Aqua; + + const PhaseUsage& pu = baseif_.phaseUsage(); + assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); + const int oil_pos = pu.phase_pos[Oil]; + + auto& segments = well_state.segments(baseif_.indexOfWell()); + auto& segment_rates = segments.rates; + auto& segment_pressure = segments.pressure; + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + std::vector fractions(baseif_.numPhases(), 0.0); + fractions[oil_pos] = 1.0; + + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int water_pos = pu.phase_pos[Water]; + fractions[water_pos] = primary_variables_[seg][WFrac]; + fractions[oil_pos] -= fractions[water_pos]; + } + + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + 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 < baseif_.numPhases(); ++p) { + const double scale = baseif_.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]; + for (int p = 0; p < baseif_.numPhases(); ++p) { + const double phase_rate = g_total * fractions[p]; + segment_rates[seg*baseif_.numPhases() + p] = phase_rate; + if (seg == 0) { // top segment + well_state.wellRates(baseif_.indexOfWell())[p] = phase_rate; + } + } + + // update the segment pressure + segment_pressure[seg] = primary_variables_[seg][SPres]; + if (seg == 0) { // top segment + well_state.update_bhp(baseif_.indexOfWell(), segment_pressure[seg]); + } + } + updateThp(well_state, rho, deferred_logger); +} + +template +void +MultisegmentWellEval:: +assembleICDPressureEq(const int seg, + const UnitSystem& unit_system, + WellState& well_state, + DeferredLogger& deferred_logger) const +{ + // TODO: upwinding needs to be taken care of + // top segment can not be a spiral ICD device + assert(seg != 0); + + // the pressure equation is something like + // p_seg - deltaP - p_outlet = 0. + // the major part is how to calculate the deltaP + + EvalWell pressure_equation = getSegmentPressure(seg); + + EvalWell icd_pressure_drop; + switch(this->segmentSet()[seg].segmentType()) { + case Segment::SegmentType::SICD : + icd_pressure_drop = pressureDropSpiralICD(seg); + break; + case Segment::SegmentType::AICD : + icd_pressure_drop = pressureDropAutoICD(seg, unit_system); + break; + case Segment::SegmentType::VALVE : + icd_pressure_drop = pressureDropValve(seg); + break; + default: { + OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(this->segmentSet()[seg].segmentNumber()) + + " for well " + baseif_.name() + " is not of ICD type", deferred_logger); + } + } + pressure_equation = pressure_equation - icd_pressure_drop; + well_state.segments(baseif_.indexOfWell()).pressure_drop_friction[seg] = icd_pressure_drop.value(); + + const int seg_upwind = upwinding_segments_[seg]; + resWell_[seg][SPres] = pressure_equation.value(); + duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); + duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + Indices::numEq); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); + } + + // contribution from the outlet segment + const int outlet_segment_index = this->segmentNumberToIndex(this->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 + Indices::numEq); + } +} + +template +void +MultisegmentWellEval:: +assemblePressureEq(const int seg, + const UnitSystem& unit_system, + WellState& well_state, + DeferredLogger& deferred_logger) const +{ + switch(this->segmentSet()[seg].segmentType()) { + case Segment::SegmentType::SICD : + case Segment::SegmentType::AICD : + case Segment::SegmentType::VALVE : { + assembleICDPressureEq(seg, unit_system, well_state,deferred_logger); + break; + } + default : + assembleDefaultPressureEq(seg, well_state); + } +} + +template +std::vector +MultisegmentWellEval:: +getWellResiduals(const std::vector& B_avg, + DeferredLogger& deferred_logger) const +{ + assert(int(B_avg.size() ) == baseif_.numComponents()); + std::vector residuals(numWellEq + 1, 0.0); + + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + double residual = 0.; + if (eq_idx < baseif_.numComponents()) { + residual = std::abs(resWell_[seg][eq_idx]) * B_avg[eq_idx]; + } else { + if (seg > 0) { + residual = std::abs(resWell_[seg][eq_idx]); + } + } + if (std::isnan(residual) || std::isinf(residual)) { + OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for residal get for well " << baseif_.name() + << " segment " << seg << " eq_idx " << eq_idx, deferred_logger); + } + + if (residual > residuals[eq_idx]) { + residuals[eq_idx] = residual; + } + } + } + + // handling the control equation residual + { + const double control_residual = std::abs(resWell_[0][numWellEq - 1]); + if (std::isnan(control_residual) || std::isinf(control_residual)) { + OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for control residal get for well " << baseif_.name(), deferred_logger); + } + residuals[numWellEq] = control_residual; + } + + return residuals; +} + +template +double +MultisegmentWellEval:: +getControlTolerance(const WellState& well_state, + const double tolerance_wells, + const double tolerance_pressure_ms_wells, + DeferredLogger& deferred_logger) const +{ + double control_tolerance = 0.; + + const int well_index = baseif_.indexOfWell(); + if (baseif_.isInjector() ) + { + auto current = well_state.currentInjectionControl(well_index); + switch(current) { + case Well::InjectorCMode::THP: + control_tolerance = tolerance_pressure_ms_wells; + break; + case Well::InjectorCMode::BHP: + control_tolerance = tolerance_wells; + break; + case Well::InjectorCMode::RATE: + case Well::InjectorCMode::RESV: + control_tolerance = tolerance_wells; + break; + case Well::InjectorCMode::GRUP: + control_tolerance = tolerance_wells; + break; + default: + OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger); + } + } + + if (baseif_.isProducer() ) + { + auto current = well_state.currentProductionControl(well_index); + switch(current) { + case Well::ProducerCMode::THP: + control_tolerance = tolerance_pressure_ms_wells; // 0.1 bar + break; + case Well::ProducerCMode::BHP: + control_tolerance = tolerance_wells; // 0.01 bar + break; + case Well::ProducerCMode::ORAT: + case Well::ProducerCMode::WRAT: + case Well::ProducerCMode::GRAT: + case Well::ProducerCMode::LRAT: + case Well::ProducerCMode::RESV: + case Well::ProducerCMode::CRAT: + control_tolerance = tolerance_wells; // smaller tolerance for rate control + break; + case Well::ProducerCMode::GRUP: + control_tolerance = tolerance_wells; // smaller tolerance for rate control + break; + default: + OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger); + } + } + + return control_tolerance; +} + +template +double +MultisegmentWellEval:: +getResidualMeasureValue(const WellState& well_state, + const std::vector& residuals, + const double tolerance_wells, + const double tolerance_pressure_ms_wells, + DeferredLogger& deferred_logger) const +{ + assert(int(residuals.size()) == numWellEq + 1); + + const double rate_tolerance = tolerance_wells; + int count = 0; + double sum = 0; + for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) { + if (residuals[eq_idx] > rate_tolerance) { + sum += residuals[eq_idx] / rate_tolerance; + ++count; + } + } + + const double pressure_tolerance = tolerance_pressure_ms_wells; + if (residuals[SPres] > pressure_tolerance) { + sum += residuals[SPres] / pressure_tolerance; + ++count; + } + + const double control_tolerance = getControlTolerance(well_state, + tolerance_wells, + tolerance_pressure_ms_wells, + deferred_logger); + if (residuals[SPres + 1] > control_tolerance) { + sum += residuals[SPres + 1] / control_tolerance; + ++count; + } + + // if (count == 0), it should be converged. + assert(count != 0); + + return sum; +} + +template +void +MultisegmentWellEval:: +updateUpwindingSegments() +{ + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + // special treatment is needed for segment 0 + if (seg == 0) { + // we are not supposed to have injecting producers and producing injectors + assert( ! (baseif_.isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) ); + assert( ! (baseif_.isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) ); + upwinding_segments_[seg] = seg; + continue; + } + + // for other normal segments + if (primary_variables_evaluation_[seg][GTotal] <= 0.) { + upwinding_segments_[seg] = seg; + } else { + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); + upwinding_segments_[seg] = outlet_segment_index; + } + } +} + +#if HAVE_CUDA || HAVE_OPENCL +template +void +MultisegmentWellEval:: +addWellContribution(WellContributions& wellContribs) const +{ + unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_ + unsigned int BnumBlocks = duneB_.nonzeroes(); + unsigned int DnumBlocks = duneD_.nonzeroes(); + + // duneC + std::vector Ccols; + std::vector Cvals; + Ccols.reserve(BnumBlocks); + Cvals.reserve(BnumBlocks * Indices::numEq * numWellEq); + for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) { + for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) { + Ccols.emplace_back(colC.index()); + for (int i = 0; i < numWellEq; ++i) { + for (int j = 0; j < Indices::numEq; ++j) { + Cvals.emplace_back((*colC)[i][j]); + } + } + } + } + + // duneD + Dune::UMFPack umfpackMatrix(duneD_, 0); + double *Dvals = umfpackMatrix.getInternalMatrix().getValues(); + auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart(); + auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex(); + + // duneB + std::vector Bcols; + std::vector Brows; + std::vector Bvals; + Bcols.reserve(BnumBlocks); + Brows.reserve(Mb+1); + Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq); + Brows.emplace_back(0); + unsigned int sumBlocks = 0; + for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) { + int sizeRow = 0; + for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) { + Bcols.emplace_back(colB.index()); + for (int i = 0; i < numWellEq; ++i) { + for (int j = 0; j < Indices::numEq; ++j) { + Bvals.emplace_back((*colB)[i][j]); + } + } + sizeRow++; + } + sumBlocks += sizeRow; + Brows.emplace_back(sumBlocks); + } + + wellContribs.addMultisegmentWellContribution(Indices::numEq, + numWellEq, + Mb, + Bvals, + Bcols, + Brows, + DnumBlocks, + Dvals, + Dcols, + Drows, + Cvals); +} +#endif + +#define INSTANCE(A,...) \ +template class MultisegmentWellEval,__VA_ARGS__,double>; + +// One phase +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u>) + +// Two phase +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>) + +// Blackoil +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u>) + +// Alt indices +INSTANCE(EclAlternativeBlackOilIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) + +} // namespace Opm diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp new file mode 100644 index 000000000..ad2e69ce3 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -0,0 +1,318 @@ +/* + 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_EVAL_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED + +#include + +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm +{ + +class ConvergenceReport; +class GroupState; +class Schedule; +class WellContributions; +template class WellInterfaceIndices; +class WellState; + +template +class MultisegmentWellEval : public MultisegmentWellGeneric +{ +public: +#if HAVE_CUDA || HAVE_OPENCL + /// add the contribution (C, D, B matrices) of this Well to the WellContributions object + void addWellContribution(WellContributions& wellContribs) const; +#endif + +protected: + // 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. + + // Table showing the primary variable indices, depending on what phases are present: + // + // WOG OG WG WO W/O/G (single phase) + // GTotal 0 0 0 0 0 + // WFrac 1 -1000 1 1 -1000 + // GFrac 2 1 -1000 -1000 -1000 + // Spres 3 2 2 2 1 + + static constexpr bool has_water = (Indices::waterSaturationIdx >= 0); + static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0); + static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0; + + // In the implementation, one should use has_wfrac_variable + // rather than has_water to check if you should do something + // with the variable at the WFrac location, similar for GFrac. + static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1; + static constexpr bool has_gfrac_variable = has_gas && has_oil; + + static constexpr int GTotal = 0; + static constexpr int WFrac = has_wfrac_variable ? 1 : -1000; + static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000; + static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1; + + // the number of well equations TODO: it should have a more general strategy for it + static constexpr int numWellEq = Indices::numPhases + 1; + + // 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 + using VectorBlockWellType = Dune::FieldVector; + using BVectorWell = Dune::BlockVector; + + using VectorBlockType = Dune::FieldVector; + using BVector = Dune::BlockVector; + + // the matrix type for the diagonal matrix D + using DiagMatrixBlockWellType = Dune::FieldMatrix; + using DiagMatWell = Dune::BCRSMatrix; + + // the matrix type for the non-diagonal matrix B and C^T + using OffDiagMatrixBlockWellType = Dune::FieldMatrix; + using OffDiagMatWell = Dune::BCRSMatrix; + + // 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. + using EvalWell = DenseAd::Evaluation; + using Eval = DenseAd::Evaluation; + + MultisegmentWellEval(WellInterfaceIndices& baseif); + + void initMatrixAndVectors(const int num_cells) const; + void initPrimaryVariablesEvaluation() const; + + void assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + const double rho, + DeferredLogger& deferred_logger); + + void assembleDefaultPressureEq(const int seg, + WellState& well_state) const; + + + // assemble pressure equation for ICD segments + void assembleICDPressureEq(const int seg, + const UnitSystem& unit_system, + WellState& well_state, + DeferredLogger& deferred_logger) const; + + + void assemblePressureEq(const int seg, + const UnitSystem& unit_system, + WellState& well_state, + DeferredLogger& deferred_logger) const; + + void checkConvergenceControlEq(const WellState& well_state, + ConvergenceReport& report, + const double tolerance_pressure_ms_wells, + const double tolerance_wells, + const double max_residual_allowed, + DeferredLogger& deferred_logger) const; + + void computePerfRatePressure(const EvalWell& pressure_cell, + const EvalWell& rs, + const EvalWell& rv, + const std::vector& b_perfcells, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const; + + /// check whether the well equations get converged for this well + ConvergenceReport getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + DeferredLogger& deferred_logger, + const double max_residual_allowed, + const double tolerance_wells, + const double relaxed_inner_tolerance_flow_ms_well, + const double tolerance_pressure_ms_wells, + const double relaxed_inner_tolerance_pressure_ms_well, + const bool relax_tolerance) const; + + // handling the overshooting and undershooting of the fractions + void processFractions(const int seg) const; + + // xw = inv(D)*(rw - C*x) + void recoverSolutionWell(const BVector& x, + BVectorWell& xw) const; + + void updatePrimaryVariables(const WellState& well_state) const; + + void updateUpwindingSegments(); + + // updating the well_state based on well solution dwells + void updateWellState(const BVectorWell& dwells, + const double relaxation_factor, + const double DFLimit, + const double max_pressure_change) const; + + void computeSegmentFluidProperties(const EvalWell& temperature, + const EvalWell& saltConcentration, + int pvt_region_index); + + EvalWell getBhp() const; + EvalWell getFrictionPressureLoss(const int seg) const; + EvalWell getHydroPressureLoss(const int seg) const; + EvalWell getQs(const int comp_idx) const; + EvalWell getSegmentGTotal(const int seg) const; + EvalWell getSegmentPressure(const int seg) const; + EvalWell getSegmentRate(const int seg, const int comp_idx) const; + EvalWell getSegmentRateUpwinding(const int seg, + const size_t comp_idx) const; + EvalWell getSegmentSurfaceVolume(const EvalWell& temperature, + const EvalWell& saltConcentration, + const int pvt_region_index, + const int seg_idx) const; + EvalWell getWQTotal() const; + + + std::vector getWellResiduals(const std::vector& B_avg, + DeferredLogger& deferred_logger) const; + + double getControlTolerance(const WellState& well_state, + const double tolerance_wells, + const double tolerance_pressure_ms_wells, + DeferredLogger& deferred_logger) const; + + double getResidualMeasureValue(const WellState& well_state, + const std::vector& residuals, + const double tolerance_wells, + const double tolerance_pressure_ms_wells, + DeferredLogger& deferred_logger) const; + + void handleAccelerationPressureLoss(const int seg, + WellState& well_state) const; + + // pressure drop for Autonomous ICD segment (WSEGAICD) + EvalWell pressureDropAutoICD(const int seg, + const UnitSystem& unit_system) const; + + // pressure drop for Spiral ICD segment (WSEGSICD) + EvalWell pressureDropSpiralICD(const int seg) const; + + // pressure drop for sub-critical valve (WSEGVALV) + EvalWell pressureDropValve(const int seg) const; + + void updateThp(WellState& well_state, + const double rho, + DeferredLogger& deferred_logger) const; + + void updateWellStateFromPrimaryVariables(WellState& well_state, + const double rho, + DeferredLogger& deferred_logger) const; + + // 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 unsigned compIdx) 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; + + // convert a Eval from reservoir to contain the derivative related to wells + EvalWell extendEval(const Eval& in) const; + + const WellInterfaceIndices& baseif_; + + // 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. It has offdiagonal entries for inlets and outlets. + mutable DiagMatWell duneD_; + + /// \brief solver for diagonal matrix + /// + /// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials... + mutable std::shared_ptr > duneDSolver_; + + // 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_; + + // the upwinding segment for each segment based on the flow direction + std::vector upwinding_segments_; + + // the densities of segment fluids + // we should not have this member variable + std::vector segment_densities_; + + // the mass rate of the segments + std::vector segment_mass_rates_; + + // the viscosity of the segments + std::vector segment_viscosities_; + + std::vector> segment_phase_densities_; + std::vector> segment_phase_fractions_; + std::vector> segment_phase_viscosities_; + + // 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_; +}; + +} + +#endif // OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 2b333e631..07e8ee93f 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -48,17 +48,8 @@ namespace Opm const int index_of_well, const std::vector& perf_data) : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data) - , MultisegmentWellGeneric(static_cast(*this)) - , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) - , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0) + , MSWEval(static_cast&>(*this)) , segment_fluid_initial_(this->numberOfSegments(), std::vector(num_components_, 0.0)) - , segment_densities_(this->numberOfSegments(), 0.0) - , segment_viscosities_(this->numberOfSegments(), 0.0) - , segment_mass_rates_(this->numberOfSegments(), 0.0) - , upwinding_segments_(this->numberOfSegments(), 0) - , segment_phase_fractions_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? - , segment_phase_viscosities_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? - , segment_phase_densities_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? { // not handling solvent or polymer for now with multisegment well if constexpr (has_solvent) { @@ -107,12 +98,12 @@ namespace Opm // \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); + this->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]; + this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf]; } } @@ -120,96 +111,31 @@ namespace Opm - 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 = this->numberOfSegments(); - for (const std::vector& inlets : this->segment_inlets_) { - nnz_d += 2 * inlets.size(); - } - duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d); - } - duneB_.setSize(this->numberOfSegments(), num_cells, number_of_perforations_); - duneC_.setSize(this->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 = this->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 = this->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 : this->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 : this->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 : this->segment_perforations_[row.index()]) { - const int cell_idx = well_cells_[perf]; - row.insert(cell_idx); - } - } - - resWell_.resize( this->numberOfSegments() ); - - primary_variables_.resize(this->numberOfSegments()); - primary_variables_evaluation_.resize(this->numberOfSegments()); - } - - - - - template void MultisegmentWell:: initPrimaryVariablesEvaluation() const { - for (int seg = 0; seg < this->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); - } - } + this->MSWEval::initPrimaryVariablesEvaluation(); } + + + template + void + MultisegmentWell:: + updatePrimaryVariables(const WellState& well_state, DeferredLogger& /* deferred_logger */) const + { + this->MSWEval::updatePrimaryVariables(well_state); + } + + + + + + template void MultisegmentWell:: @@ -231,75 +157,20 @@ namespace Opm template ConvergenceReport MultisegmentWell:: - getWellConvergence(const WellState& well_state, const std::vector& B_avg, DeferredLogger& deferred_logger, const bool relax_tolerance) const + getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + DeferredLogger& deferred_logger, + const bool relax_tolerance) const { - assert(int(B_avg.size()) == num_components_); - - // checking if any residual is NaN or too large. The two large one is only handled for the well flux - std::vector> abs_residual(this->numberOfSegments(), std::vector(numWellEq, 0.0)); - for (int seg = 0; seg < this->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 eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - if (eq_idx < num_components_) { // phase or component mass equations - const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx]; - if (flux_residual > maximum_residual[eq_idx]) { - maximum_residual[eq_idx] = flux_residual; - } - } else { // pressure or control equation - // for the top segment (seg == 0), it is control equation, will be checked later separately - if (seg > 0) { - // Pressure equation - const double pressure_residual = abs_residual[seg][eq_idx]; - if (pressure_residual > maximum_residual[eq_idx]) { - maximum_residual[eq_idx] = pressure_residual; - } - } - } - } - } - - using CR = ConvergenceReport; - for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - if (eq_idx < num_components_) { // phase or component mass equations - const double flux_residual = maximum_residual[eq_idx]; - // TODO: the report can not handle the segment number yet. - - if (std::isnan(flux_residual)) { - report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, name()}); - } else if (flux_residual > param_.max_residual_allowed_) { - report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, name()}); - } else if (!relax_tolerance && flux_residual > param_.tolerance_wells_) { - report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()}); - } else if (flux_residual > param_.relaxed_inner_tolerance_flow_ms_well_) { - report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()}); - } - } else { // pressure equation - const double pressure_residual = maximum_residual[eq_idx]; - const int dummy_component = -1; - if (std::isnan(pressure_residual)) { - report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()}); - } else if (std::isinf(pressure_residual)) { - report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, name()}); - } else if (!relax_tolerance && pressure_residual > param_.tolerance_pressure_ms_wells_) { - report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()}); - } else if (pressure_residual > param_.relaxed_inner_tolerance_pressure_ms_well_) { - report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()}); - } - } - } - - checkConvergenceControlEq(well_state, report, deferred_logger); - - return report; + return this->MSWEval::getWellConvergence(well_state, + B_avg, + deferred_logger, + param_.max_residual_allowed_, + param_.tolerance_wells_, + param_.relaxed_inner_tolerance_flow_ms_well_, + param_.tolerance_pressure_ms_wells_, + param_.relaxed_inner_tolerance_pressure_ms_well_, + relax_tolerance); } @@ -318,15 +189,15 @@ namespace Opm // Contributions are already in the matrix itself return; } - BVectorWell Bx(duneB_.N()); + BVectorWell Bx(this->duneB_.N()); - duneB_.mv(x, Bx); + this->duneB_.mv(x, Bx); // invDBx = duneD^-1 * Bx_ - const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx); + const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx); // Ax = Ax - duneC_^T * invDBx - duneC_.mmtv(invDBx,Ax); + this->duneC_.mmtv(invDBx,Ax); } @@ -341,74 +212,13 @@ namespace Opm if (!this->isOperable() && !this->wellIsStopped()) return; // invDrw_ = duneD^-1 * resWell_ - const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); + const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_); // r = r - duneC_^T * invDrw - duneC_.mmtv(invDrw, r); + this->duneC_.mmtv(invDrw, r); } -#if HAVE_CUDA || HAVE_OPENCL - template - void - MultisegmentWell:: - addWellContribution(WellContributions& wellContribs) const - { - unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_ - unsigned int BnumBlocks = duneB_.nonzeroes(); - unsigned int DnumBlocks = duneD_.nonzeroes(); - - // duneC - std::vector Ccols; - std::vector Cvals; - Ccols.reserve(BnumBlocks); - Cvals.reserve(BnumBlocks * numEq * numWellEq); - for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) { - for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) { - Ccols.emplace_back(colC.index()); - for (int i = 0; i < numWellEq; ++i) { - for (int j = 0; j < numEq; ++j) { - Cvals.emplace_back((*colC)[i][j]); - } - } - } - } - - // duneD - Dune::UMFPack umfpackMatrix(duneD_, 0); - double *Dvals = umfpackMatrix.getInternalMatrix().getValues(); - auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart(); - auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex(); - - // duneB - std::vector Bcols; - std::vector Brows; - std::vector Bvals; - Bcols.reserve(BnumBlocks); - Brows.reserve(Mb+1); - Bvals.reserve(BnumBlocks * numEq * numWellEq); - Brows.emplace_back(0); - unsigned int sumBlocks = 0; - for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) { - int sizeRow = 0; - for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) { - Bcols.emplace_back(colB.index()); - for (int i = 0; i < numWellEq; ++i) { - for (int j = 0; j < numEq; ++j) { - Bvals.emplace_back((*colB)[i][j]); - } - } - sizeRow++; - } - sumBlocks += sizeRow; - Brows.emplace_back(sumBlocks); - } - - wellContribs.addMultisegmentWellContribution(numEq, numWellEq, Mb, Bvals, Bcols, Brows, DnumBlocks, Dvals, Dcols, Drows, Cvals); - } -#endif - - template void MultisegmentWell:: @@ -419,7 +229,7 @@ namespace Opm if (!this->isOperable() && !this->wellIsStopped()) return; BVectorWell xw(1); - recoverSolutionWell(x, xw); + this->recoverSolutionWell(x, xw); updateWellState(xw, well_state, deferred_logger); } @@ -611,100 +421,6 @@ namespace Opm - template - void - MultisegmentWell:: - updatePrimaryVariables(const WellState& well_state, DeferredLogger& /* deferred_logger */) const - { - // TODO: to test using rate conversion coefficients to see if it will be better than - // this default one - if (!this->isOperable() && !this->wellIsStopped()) return; - - const Well& well = Base::wellEcl(); - - // the index of the top segment in the WellState - const auto& segments = well_state.segments(this->index_of_well_); - const auto& segment_rates = segments.rates; - const auto& segment_pressure = segments.pressure; - const PhaseUsage& pu = phaseUsage(); - - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - // calculate the total rate for each segment - double total_seg_rate = 0.0; - // the segment pressure - primary_variables_[seg][SPres] = segment_pressure[seg]; - // 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 + p]; - } - - primary_variables_[seg][GTotal] = total_seg_rate; - if (std::abs(total_seg_rate) > 0.) { - if (has_wfrac_variable) { - const int water_pos = pu.phase_pos[Water]; - primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg + water_pos] / total_seg_rate; - } - if (has_gfrac_variable) { - const int gas_pos = pu.phase_pos[Gas]; - primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg + gas_pos] / total_seg_rate; - } - } else { // total_seg_rate == 0 - if (this->isInjector()) { - // only single phase injection handled - auto phase = well.getInjectionProperties().injectorType; - - if (has_wfrac_variable) { - if (phase == InjectorType::WATER) { - primary_variables_[seg][WFrac] = 1.0; - } else { - primary_variables_[seg][WFrac] = 0.0; - } - } - - if (has_gfrac_variable) { - if (phase == InjectorType::GAS) { - primary_variables_[seg][GFrac] = 1.0; - } else { - primary_variables_[seg][GFrac] = 0.0; - } - } - - } else if (this->isProducer()) { // producers - if (has_wfrac_variable) { - primary_variables_[seg][WFrac] = 1.0 / number_of_phases_; - } - - if (has_gfrac_variable) { - primary_variables_[seg][GFrac] = 1.0 / number_of_phases_; - } - } - } - } - } - - - - - - template - void - MultisegmentWell:: - recoverSolutionWell(const BVector& x, BVectorWell& xw) const - { - if (!this->isOperable() && !this->wellIsStopped()) return; - - BVectorWell resWell = resWell_; - // resWell = resWell - B * x - duneB_.mmv(x, resWell); - // xw = D^-1 * resWell - xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); - } - - - - - template void MultisegmentWell:: @@ -714,7 +430,7 @@ namespace Opm // 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::applyUMFPack(duneD_, duneDSolver_, resWell_); + const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_); updateWellState(dx_well, well_state, deferred_logger); } @@ -770,7 +486,7 @@ namespace Opm } average_density /= sum_kr; - cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf]; + this->cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * this->cell_perforation_depth_diffs_[perf]; } } @@ -787,7 +503,7 @@ namespace Opm // TODO: trying to reduce the times for the surfaceVolumeFraction calculation const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value(); for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - segment_fluid_initial_[seg][comp_idx] = surface_volume * surfaceVolumeFraction(seg, comp_idx).value(); + segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value(); } } } @@ -808,48 +524,12 @@ namespace Opm 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_; + this->MSWEval::updateWellState(dwells, + relaxation_factor, + dFLimit, + max_pressure_change); - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - if (has_wfrac_variable) { - 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 (has_gfrac_variable) { - 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] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); - } - - // 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]; - - // make sure that no injector produce and no producer inject - if (seg == 0) { - if (this->isInjector()) { - primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0); - } else { - primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0); - } - } - } - - } - - updateWellStateFromPrimaryVariables(well_state, deferred_logger); + this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger); Base::calculateReservoirRates(well_state); } @@ -949,7 +629,7 @@ namespace Opm MultisegmentWell:: addWellContributions(SparseMatrixAdapter& jacobian) const { - const auto invDuneD = mswellhelpers::invertWithUMFPack(duneD_, duneDSolver_); + const auto invDuneD = mswellhelpers::invertWithUMFPack(this->duneD_, this->duneDSolver_); // We need to change matrix A as follows // A -= C^T D^-1 B @@ -959,11 +639,11 @@ namespace Opm // perforation at cell j connected to segment i. The code // assumes that no cell is connected to more than one segment, // i.e. the columns of B/C have no more than one nonzero. - for (size_t rowC = 0; rowC < duneC_.N(); ++rowC) { - for (auto colC = duneC_[rowC].begin(), endC = duneC_[rowC].end(); colC != endC; ++colC) { + for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) { + for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) { const auto row_index = colC.index(); - for (size_t rowB = 0; rowB < duneB_.N(); ++rowB) { - for (auto colB = duneB_[rowB].begin(), endB = duneB_[rowB].end(); colB != endB; ++colB) { + for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) { + for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) { const auto col_index = colB.index(); OffDiagMatrixBlockWellType tmp1; Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); @@ -980,77 +660,6 @@ namespace Opm - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - volumeFraction(const int seg, const unsigned compIdx) const - { - - if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { - return primary_variables_evaluation_[seg][WFrac]; - } - - if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { - return primary_variables_evaluation_[seg][GFrac]; - } - - // Oil fraction - EvalWell oil_fraction = 1.0; - if (has_wfrac_variable) { - oil_fraction -= primary_variables_evaluation_[seg][WFrac]; - } - - if (has_gfrac_variable) { - 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(ebosCompIdxToFlowCompIdx(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.; - for (int idx = 0; idx < num_components_; ++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:: @@ -1068,18 +677,11 @@ namespace Opm DeferredLogger& deferred_logger) const { - std::vector cmix_s(num_components_, 0.0); - - // the composition of the components inside wellbore - for (int comp_idx = 0; comp_idx < num_components_; ++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()); + const EvalWell pressure_cell = this->extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + const EvalWell rs = this->extendEval(fs.Rs()); + const EvalWell rv = this->extendEval(fs.Rv()); // not using number_of_phases_ because of solvent std::vector b_perfcells(num_components_, 0.0); @@ -1090,137 +692,24 @@ namespace Opm } const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - b_perfcells[compIdx] = extendEval(fs.invB(phaseIdx)); + b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx)); } - // pressure difference between the segment and the perforation - const EvalWell perf_seg_press_diff = gravity_ * this->segment_densities_[seg] * this->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]; - - perf_press = pressure_cell - cell_perf_press_diff; - // 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 = perf_press - (segment_pressure + perf_seg_press_diff); - - // producing perforations - if ( drawdown > 0.0) { - // Do nothing is crossflow is not allowed - if (!allow_cf && this->isInjector()) { - return; - } - - // compute component volumetric rates at standard conditions - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); - cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const EvalWell cq_s_oil = cq_s[oilCompIdx]; - const EvalWell cq_s_gas = cq_s[gasCompIdx]; - cq_s[gasCompIdx] += rs * cq_s_oil; - cq_s[oilCompIdx] += rv * cq_s_gas; - } - } else { // injecting perforations - // Do nothing if crossflow is not allowed - if (!allow_cf && this->isProducer()) { - return; - } - - // for injecting perforations, we use total mobility - EvalWell total_mob = mob_perfcells[0]; - for (int comp_idx = 1; comp_idx < num_components_; ++comp_idx) { - total_mob += mob_perfcells[comp_idx]; - } - - // injection perforations total volume rates - const EvalWell cqt_i = - Tw * (total_mob * drawdown); - - // compute volume ratio between connection and at standard conditions - EvalWell volume_ratio = 0.0; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // 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_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << name() << " during flux calcuation" - << " with rs " << rs << " and rv " << rv, deferred_logger); - } - - const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; - volume_ratio += tmp_oil / b_perfcells[oilCompIdx]; - - const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d; - volume_ratio += tmp_gas / b_perfcells[gasCompIdx]; - } else { // not having gas and oil at the same time - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx]; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx]; - } - } - // injecting connections total volumerates at standard conditions - EvalWell cqt_is = cqt_i / volume_ratio; - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; - } - } // end for injection perforations - - // calculating the perforation solution gas rate and solution oil rates - if (this->isProducer()) { - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells - // s means standard condition, r means reservoir condition - // q_os = q_or * b_o + rv * q_gr * b_g - // q_gs = q_gr * g_g + rs * q_or * b_o - // d = 1.0 - rs * rv - // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) - // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) - - const double d = 1.0 - rv.value() * rs.value(); - // vaporized oil into gas - // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d - perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d; - // dissolved of gas in oil - // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d - perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d; - } - } - } - - - - - - 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; + this->MSWEval::computePerfRatePressure(pressure_cell, + rs, + rv, + b_perfcells, + mob_perfcells, + Tw, + seg, + perf, + segment_pressure, + allow_cf, + cq_s, + perf_press, + perf_dis_gas_rate, + perf_vap_oil_rate, + deferred_logger); } @@ -1253,279 +742,13 @@ namespace Opm const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value()); - saltConcentration = extendEval(fs.saltConcentration()); + saltConcentration = this->extendEval(fs.saltConcentration()); pvt_region_index = fs.pvtRegionIndex(); } - std::vector surf_dens(num_components_); - // Surface density. - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index ); - } - - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - // the compostion of the components inside wellbore under surface condition - std::vector mix_s(num_components_, 0.0); - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); - } - - std::vector b(num_components_, 0.0); - std::vector visc(num_components_, 0.0); - std::vector& phase_densities = this->segment_phase_densities_[seg]; - - const EvalWell seg_pressure = getSegmentPressure(seg); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - b[waterCompIdx] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration); - visc[waterCompIdx] = - FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration); - // TODO: double check here - // TODO: should not we use phaseIndex here? - phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx]; - } - - EvalWell rv(0.0); - // gas phase - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); - if (mix_s[oilCompIdx] > 0.0) { - if (mix_s[gasCompIdx] > 0.0) { - rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; - } - - if (rv > rvmax) { - rv = rvmax; - } - b[gasCompIdx] = - FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); - visc[gasCompIdx] = - FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv); - phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx] - + rv * b[gasCompIdx] * surf_dens[oilCompIdx]; - } else { // no oil exists - b[gasCompIdx] = - FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[gasCompIdx] = - FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); - phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx]; - } - } else { // no Liquid phase - // it is the same with zero mix_s[Oil] - b[gasCompIdx] = - FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[gasCompIdx] = - FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); - } - } - - EvalWell rs(0.0); - // oil phase - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); - if (mix_s[gasCompIdx] > 0.0) { - if (mix_s[oilCompIdx] > 0.0) { - rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; - } - - if (rs > rsmax) { - rs = rsmax; - } - b[oilCompIdx] = - FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); - visc[oilCompIdx] = - FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs); - phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx] - + rs * b[oilCompIdx] * surf_dens[gasCompIdx]; - } else { // no oil exists - b[oilCompIdx] = - FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[oilCompIdx] = - FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); - phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx]; - } - } else { // no Liquid phase - // it is the same with zero mix_s[Oil] - b[oilCompIdx] = - FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - visc[oilCompIdx] = - FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure); - } - } - - segment_phase_viscosities_[seg] = visc; - - std::vector mix(mix_s); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - - const EvalWell d = 1.0 - rs * rv; - - if (rs != 0.0) { // rs > 0.0? - mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; - } - if (rv != 0.0) { // rv > 0.0? - mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; - } - } - - EvalWell volrat(0.0); - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - volrat += mix[comp_idx] / b[comp_idx]; - } - - segment_viscosities_[seg] = 0.; - // calculate the average viscosity - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; - // TODO: a little more work needs to be done to handle the negative fractions here - segment_phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; - segment_viscosities_[seg] += visc[comp_idx] * segment_phase_fractions_[seg][comp_idx]; - } - - EvalWell density(0.0); - for (int comp_idx = 0; comp_idx < num_components_; ++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 comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell rate = getSegmentRateUpwinding(seg, comp_idx); - segment_mass_rates_[seg] += rate * surf_dens[comp_idx]; - } - } - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getSegmentPressure(const int seg) const - { - return primary_variables_evaluation_[seg][SPres]; - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getBhp() const - { - return getSegmentPressure(0); - } - - - - - - 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:: - getQs(const int comp_idx) const - { - return getSegmentRate(0, comp_idx); - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getSegmentRateUpwinding(const int seg, - const size_t comp_idx) const - { - const int seg_upwind = upwinding_segments_[seg]; - // the result will contain the derivative with resepct to GTotal in segment seg, - // and the derivatives with respect to WFrac GFrac in segment seg_upwind. - // the derivative with respect to SPres should be zero. - if (seg == 0 && this->isInjector()) { - const Well& well = Base::wellEcl(); - auto phase = well.getInjectionProperties().injectorType; - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) - && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx - && phase == InjectorType::WATER) - return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx)); - - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) - && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx - && phase == InjectorType::OIL) - return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx)); - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) - && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx - && phase == InjectorType::GAS) - return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx)); - - return 0.0; - } - - const EvalWell segment_rate = primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx); - - assert(segment_rate.derivative(SPres + numEq) == 0.); - - return segment_rate; - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getSegmentGTotal(const int seg) const - { - return primary_variables_evaluation_[seg][GTotal]; - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getWQTotal() const - { - return getSegmentGTotal(0); + this->MSWEval::computeSegmentFluidProperties(temperature, + saltConcentration, + pvt_region_index); } @@ -1557,7 +780,7 @@ namespace Opm } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx)); + mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx)); } // if (has_solvent) { // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); @@ -1578,358 +801,22 @@ namespace Opm } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); + mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); } - } } - template - void - MultisegmentWell:: - assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::InjectionControls& inj_controls, - const Well::ProductionControls& prod_controls, - DeferredLogger& deferred_logger) - { - - EvalWell control_eq(0.0); - - const auto& well = well_ecl_; - - auto getRates = [&]() { - std::vector rates(3, 0.0); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)); - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); - } - return rates; - }; - - if (this->wellIsStopped()) { - control_eq = getWQTotal(); - } else if (this->isInjector() ) { - // Find scaling factor to get injection rate, - const InjectorType injectorType = inj_controls.injector_type; - double scaling = 1.0; - const auto& pu = phaseUsage(); - switch (injectorType) { - case InjectorType::WATER: - { - scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]); - break; - } - case InjectorType::OIL: - { - scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]); - break; - } - case InjectorType::GAS: - { - scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]); - break; - } - default: - throw("Expected WATER, OIL or GAS as type for injectors " + well.name()); - } - const EvalWell injection_rate = getWQTotal() / scaling; - // Setup function for evaluation of BHP from THP (used only if needed). - auto bhp_from_thp = [&]() { - const auto rates = getRates(); - return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger); - }; - // Call generic implementation. - Base::assembleControlEqInj(well_state, group_state, schedule, summaryState, inj_controls, getBhp(), injection_rate, bhp_from_thp, control_eq, deferred_logger); - } else { - // Find rates. - const auto rates = getRates(); - // Setup function for evaluation of BHP from THP (used only if needed). - auto bhp_from_thp = [&]() { - return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger); - }; - // Call generic implementation. - Base::assembleControlEqProd(well_state, group_state, schedule, summaryState, prod_controls, getBhp(), rates, bhp_from_thp, control_eq, deferred_logger); - } - - // 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:: - updateThp(WellState& well_state, DeferredLogger& deferred_logger) const - { - // When there is no vaild VFP table provided, we set the thp to be zero. - if (!this->isVFPActive(deferred_logger) || this->wellIsStopped()) { - well_state.update_thp(index_of_well_, 0.); - return; - } - - // the well is under other control types, we calculate the thp based on bhp and rates - std::vector rates(3, 0.0); - - const PhaseUsage& pu = phaseUsage(); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - rates[ Water ] = well_state.wellRates(index_of_well_)[pu.phase_pos[ Water ] ]; - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - rates[ Oil ] = well_state.wellRates(index_of_well_)[pu.phase_pos[ Oil ] ]; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - rates[ Gas ] = well_state.wellRates(index_of_well_)[pu.phase_pos[ Gas ] ]; - } - - const double bhp = well_state.bhp(index_of_well_); - - well_state.update_thp(index_of_well_, this->calculateThpFromBhp(rates, bhp, getRefDensity(), deferred_logger)); - - } - - - template double MultisegmentWell:: getRefDensity() const { - return segment_densities_[0].value(); + return this->segment_densities_[0].value(); } - template - void - MultisegmentWell:: - assembleDefaultPressureEq(const int seg, WellState& well_state) 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 - // TODO: we might be able to add member variables to store these values, then we update well state - // after converged - const auto hydro_pressure_drop = getHydroPressureLoss(seg); - auto& segments = well_state.segments(this->index_of_well_); - segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value(); - pressure_equation -= hydro_pressure_drop; - - if (this->frictionalPressureLossConsidered()) { - const auto friction_pressure_drop = getFrictionPressureLoss(seg); - pressure_equation -= friction_pressure_drop; - segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); - } - - resWell_[seg][SPres] = pressure_equation.value(); - const int seg_upwind = upwinding_segments_[seg]; - duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq); - duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq); - if (has_wfrac_variable) { - duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq); - } - if (has_gfrac_variable) { - duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq); - } - - // contribution from the outlet segment - const int outlet_segment_index = this->segmentNumberToIndex(this->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 (this->accelerationalPressureLossConsidered()) { - handleAccelerationPressureLoss(seg, well_state); - } - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getHydroPressureLoss(const int seg) const - { - return segment_densities_[seg] * gravity_ * this->segment_depth_diffs_[seg]; - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getFrictionPressureLoss(const int seg) const - { - const EvalWell mass_rate = segment_mass_rates_[seg]; - const int seg_upwind = upwinding_segments_[seg]; - EvalWell density = segment_densities_[seg_upwind]; - EvalWell visc = segment_viscosities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (seg != seg_upwind) { - density.clearDerivatives(); - visc.clearDerivatives(); - } - const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); - const double length = this->segmentSet()[seg].totalLength() - this->segmentSet()[outlet_segment_index].totalLength(); - assert(length > 0.); - const double roughness = this->segmentSet()[seg].roughness(); - const double area = this->segmentSet()[seg].crossArea(); - const double diameter = this->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, WellState& well_state) const - { - const double area = this->segmentSet()[seg].crossArea(); - const EvalWell mass_rate = segment_mass_rates_[seg]; - const int seg_upwind = upwinding_segments_[seg]; - EvalWell density = segment_densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (seg != seg_upwind) { - density.clearDerivatives(); - } - - EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density); - // handling the velocity head of intlet segments - for (const int inlet : this->segment_inlets_[seg]) { - const int seg_upwind_inlet = upwinding_segments_[inlet]; - const double inlet_area = this->segmentSet()[inlet].crossArea(); - EvalWell inlet_density = this->segment_densities_[seg_upwind_inlet]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (inlet != seg_upwind_inlet) { - inlet_density.clearDerivatives(); - } - const EvalWell inlet_mass_rate = segment_mass_rates_[inlet]; - accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density); - } - - // We change the sign of the accelerationPressureLoss for injectors. - // Is this correct? Testing indicates that this is what the reference simulator does - const double sign = mass_rate < 0. ? 1.0 : - 1.0; - accelerationPressureLoss *= sign; - - well_state.segments(this->index_of_well_).pressure_drop_accel[seg] = accelerationPressureLoss.value(); - - resWell_[seg][SPres] -= accelerationPressureLoss.value(); - duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + numEq); - duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq); - if (has_wfrac_variable) { - duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq); - } - if (has_gfrac_variable) { - duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq); - } - } - - - - - - template - void - MultisegmentWell:: - processFractions(const int seg) const - { - const PhaseUsage& pu = phaseUsage(); - - std::vector fractions(number_of_phases_, 0.0); - - assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); - const int oil_pos = pu.phase_pos[Oil]; - fractions[oil_pos] = 1.0; - - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; - fractions[water_pos] = primary_variables_[seg][WFrac]; - fractions[oil_pos] -= fractions[water_pos]; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; - fractions[gas_pos] = primary_variables_[seg][GFrac]; - fractions[oil_pos] -= fractions[gas_pos]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int water_pos = pu.phase_pos[Water]; - if (fractions[water_pos] < 0.0) { - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]); - } - fractions[oil_pos] /= (1.0 - fractions[water_pos]); - fractions[water_pos] = 0.0; - } - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int gas_pos = pu.phase_pos[Gas]; - if (fractions[gas_pos] < 0.0) { - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - 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 ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]); - } - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]); - } - fractions[oil_pos] = 0.0; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]]; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]]; - } - } - - template void MultisegmentWell:: @@ -2006,7 +893,7 @@ namespace Opm for (int seg = 0; seg < nseg; ++seg) { // calculating the perforation rate for each perforation that belongs to this segment const double segment_depth = this->segmentSet()[seg].depth(); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, segment_densities_[seg].value(), gravity_); + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, this->segment_densities_[seg].value(), gravity_); ref_depth = segment_depth; seg_bhp_press_diff += dp; for (const int perf : this->segment_perforations_[seg]) { @@ -2021,9 +908,9 @@ namespace Opm const auto& fs = int_quantities.fluidState(); // the pressure of the reservoir grid block the well connection is in // pressure difference between the segment and the perforation - const double perf_seg_press_diff = gravity_ * segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf]; + const double perf_seg_press_diff = gravity_ * this->segment_densities_[seg].value() * this->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]; + const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value(); // calculating the b for the connection @@ -2131,68 +1018,6 @@ namespace Opm - template - void - MultisegmentWell:: - updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const - { - const PhaseUsage& pu = phaseUsage(); - assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); - const int oil_pos = pu.phase_pos[Oil]; - - auto& segments = well_state.segments(this->index_of_well_); - auto& segment_rates = segments.rates; - auto& segment_pressure = segments.pressure; - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - std::vector fractions(number_of_phases_, 0.0); - fractions[oil_pos] = 1.0; - - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; - fractions[water_pos] = primary_variables_[seg][WFrac]; - fractions[oil_pos] -= fractions[water_pos]; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - 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]; - for (int p = 0; p < number_of_phases_; ++p) { - const double phase_rate = g_total * fractions[p]; - segment_rates[seg*this->number_of_phases_ + p] = phase_rate; - if (seg == 0) { // top segment - well_state.wellRates(index_of_well_)[p] = phase_rate; - } - } - - // update the segment pressure - segment_pressure[seg] = primary_variables_[seg][SPres]; - if (seg == 0) { // top segment - well_state.update_bhp(index_of_well_, segment_pressure[seg]); - } - } - updateThp(well_state, deferred_logger); - } - - - - template bool MultisegmentWell:: @@ -2208,7 +1033,7 @@ namespace Opm const int max_iter_number = param_.max_inner_iter_ms_wells_; const WellState well_state0 = well_state; - const std::vector residuals0 = getWellResiduals(Base::B_avg_, deferred_logger); + const std::vector residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger); std::vector > residual_history; std::vector measure_history; int it = 0; @@ -2222,7 +1047,7 @@ namespace Opm assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); - const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); + const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_); if (it > param_.strict_inner_iter_ms_wells_) relax_convergence = true; @@ -2233,8 +1058,12 @@ namespace Opm break; } - residual_history.push_back(getWellResiduals(Base::B_avg_, deferred_logger)); - measure_history.push_back(getResidualMeasureValue(well_state, residual_history[it], deferred_logger) ); + residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger)); + measure_history.push_back(this->getResidualMeasureValue(well_state, + residual_history[it], + param_.tolerance_wells_, + param_.tolerance_pressure_ms_wells_, + deferred_logger) ); bool is_oscillate = false; bool is_stagnate = false; @@ -2327,19 +1156,19 @@ namespace Opm if (!this->isOperable() && !this->wellIsStopped()) return; // update the upwinding segments - updateUpwindingSegments(); + this->updateUpwindingSegments(); // calculate the fluid properties needed. computeSegmentFluidProperties(ebosSimulator); // clear all entries - duneB_ = 0.0; - duneC_ = 0.0; + this->duneB_ = 0.0; + this->duneC_ = 0.0; - duneD_ = 0.0; - resWell_ = 0.0; + this->duneD_ = 0.0; + this->resWell_ = 0.0; - duneDSolver_.reset(); + this->duneDSolver_.reset(); well_state.wellVaporizedOilRates(index_of_well_) = 0.; well_state.wellDissolvedGasRates(index_of_well_) = 0.; @@ -2365,30 +1194,30 @@ namespace Opm const Scalar regularization_factor = param_.regularization_factor_ms_wells_; // for each component for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx) + const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx) - segment_fluid_initial_[seg][comp_idx]) / dt; - resWell_[seg][comp_idx] += accumulation_term.value(); + this->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); + this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq); } } } // considering the contributions due to flowing out from the segment { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell segment_rate = getSegmentRateUpwinding(seg, comp_idx) * well_efficiency_factor_; + const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * well_efficiency_factor_; - const int seg_upwind = upwinding_segments_[seg]; + const int seg_upwind = this->upwinding_segments_[seg]; // segment_rate contains the derivatives with respect to GTotal in seg, // and WFrac and GFrac in seg_upwind - resWell_[seg][comp_idx] -= segment_rate.value(); - duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq); + this->resWell_[seg][comp_idx] -= segment_rate.value(); + this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq); + this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq); + this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq); } // pressure derivative should be zero } @@ -2398,18 +1227,18 @@ namespace Opm { for (const int inlet : this->segment_inlets_[seg]) { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell inlet_rate = getSegmentRateUpwinding(inlet, comp_idx) * well_efficiency_factor_; + const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * well_efficiency_factor_; - const int inlet_upwind = upwinding_segments_[inlet]; + const int inlet_upwind = this->upwinding_segments_[inlet]; // inlet_rate contains the derivatives with respect to GTotal in inlet, // and WFrac and GFrac in inlet_upwind - resWell_[seg][comp_idx] += inlet_rate.value(); - duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq); + this->resWell_[seg][comp_idx] += inlet_rate.value(); + this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq); + this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq); + this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq); } // pressure derivative should be zero } @@ -2417,7 +1246,7 @@ namespace Opm } // calculating the perforation rate for each perforation that belongs to this segment - const EvalWell seg_pressure = getSegmentPressure(seg); + const EvalWell seg_pressure = this->getSegmentPressure(seg); auto& perf_rates = well_state.perfPhaseRates(this->index_of_well_); auto& perf_press_state = well_state.perfPress(this->index_of_well_); for (const int perf : this->segment_perforations_[seg]) { @@ -2452,21 +1281,21 @@ namespace Opm connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); // subtract sum of phase fluxes in the well equations. - resWell_[seg][comp_idx] += cq_s_effective.value(); + this->resWell_[seg][comp_idx] += cq_s_effective.value(); // assemble the jacobians for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { // also need to consider the efficiency factor when manipulating the jacobians. - duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix + this->duneC_[seg][cell_idx][pv_idx][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); + this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + numEq); } for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { // also need to consider the efficiency factor when manipulating the jacobians. - duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); + this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); } } } @@ -2475,10 +1304,17 @@ namespace Opm if (seg == 0) { // top segment, pressure equation is the control equation const auto& summaryState = ebosSimulator.vanguard().summaryState(); const Schedule& schedule = ebosSimulator.vanguard().schedule(); - assembleControlEq(well_state, group_state, schedule, summaryState, inj_controls, prod_controls, deferred_logger); + this->assembleControlEq(well_state, + group_state, + schedule, + summaryState, + inj_controls, + prod_controls, + getRefDensity(), + deferred_logger); } else { const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem(); - assemblePressureEq(seg, unit_system, well_state, deferred_logger); + this->assemblePressureEq(seg, unit_system, well_state, deferred_logger); } } } @@ -2486,26 +1322,6 @@ namespace Opm - template - void - MultisegmentWell:: - assemblePressureEq(const int seg, const UnitSystem& unit_system, - WellState& well_state, DeferredLogger& deferred_logger) const - { - switch(this->segmentSet()[seg].segmentType()) { - case Segment::SegmentType::SICD : - case Segment::SegmentType::AICD : - case Segment::SegmentType::VALVE : { - assembleICDPressureEq(seg, unit_system, well_state,deferred_logger); - break; - } - default : - assembleDefaultPressureEq(seg, well_state); - } - } - - - template bool MultisegmentWell:: @@ -2524,7 +1340,7 @@ namespace Opm const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { - const EvalWell segment_pressure = getSegmentPressure(seg); + const EvalWell segment_pressure = this->getSegmentPressure(seg); for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; @@ -2532,9 +1348,9 @@ namespace Opm const auto& fs = intQuants.fluidState(); // pressure difference between the segment and the perforation - const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf]; + const EvalWell perf_seg_press_diff = gravity_ * this->segment_densities_[seg] * this->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]; + const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value(); const double perf_press = pressure_cell - cell_perf_press_diff; @@ -2585,450 +1401,20 @@ namespace Opm const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value()); - saltConcentration = extendEval(fs.saltConcentration()); + saltConcentration = this->extendEval(fs.saltConcentration()); pvt_region_index = fs.pvtRegionIndex(); } - const EvalWell seg_pressure = getSegmentPressure(seg_idx); - - std::vector mix_s(num_components_, 0.0); - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - mix_s[comp_idx] = surfaceVolumeFraction(seg_idx, comp_idx); - } - - std::vector b(num_components_, 0.); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - b[waterCompIdx] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration); - } - - EvalWell rv(0.0); - // gas phase - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); - if (rvmax < 0.0) { // negative rvmax can happen if the seg_pressure is outside the range of the table - rvmax = 0.0; - } - if (mix_s[oilCompIdx] > 0.0) { - if (mix_s[gasCompIdx] > 0.0) { - rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; - } - - if (rv > rvmax) { - rv = rvmax; - } - b[gasCompIdx] = - FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); - } else { // no oil exists - b[gasCompIdx] = - FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - } - } else { // no Liquid phase - // it is the same with zero mix_s[Oil] - b[gasCompIdx] = - FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - } - } - - EvalWell rs(0.0); - // oil phase - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); - if (rsmax < 0.0) { // negative rsmax can happen if the seg_pressure is outside the range of the table - rsmax = 0.0; - } - if (mix_s[gasCompIdx] > 0.0) { - if (mix_s[oilCompIdx] > 0.0) { - rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; - } - // std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl; - - if (rs > rsmax) { - rs = rsmax; - } - b[oilCompIdx] = - FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); - } else { // no oil exists - b[oilCompIdx] = - FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - } - } else { // no gas phase - // it is the same with zero mix_s[Gas] - b[oilCompIdx] = - FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); - } - } - - std::vector mix(mix_s); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - - const EvalWell d = 1.0 - rs * rv; - if (d <= 0.0 || d > 1.0) { - std::ostringstream sstr; - sstr << "Problematic d value " << d << " obtained for well " << name() - << " during conversion to surface volume with rs " << rs - << ", rv " << rv << " and pressure " << seg_pressure - << " obtaining d " << d; - OpmLog::debug(sstr.str()); - OPM_THROW_NOLOG(NumericalIssue, sstr.str()); - } - - if (rs > 0.0) { // rs > 0.0? - mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; - } - if (rv > 0.0) { // rv > 0.0? - mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; - } - } - - EvalWell vol_ratio(0.0); - for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - vol_ratio += mix[comp_idx] / b[comp_idx]; - } - - // We increase the segment volume with a factor 10 to stabilize the system. - const double volume = this->segmentSet()[seg_idx].volume(); - - return volume / vol_ratio; + return this->MSWEval::getSegmentSurfaceVolume(temperature, + saltConcentration, + pvt_region_index, + seg_idx); } - template - std::vector::Scalar> - MultisegmentWell:: - getWellResiduals(const std::vector& B_avg, - DeferredLogger& deferred_logger) const - { - assert(int(B_avg.size() ) == num_components_); - std::vector residuals(numWellEq + 1, 0.0); - - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - double residual = 0.; - if (eq_idx < num_components_) { - residual = std::abs(resWell_[seg][eq_idx]) * B_avg[eq_idx]; - } else { - if (seg > 0) { - residual = std::abs(resWell_[seg][eq_idx]); - } - } - if (std::isnan(residual) || std::isinf(residual)) { - OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for residal get for well " << name() - << " segment " << seg << " eq_idx " << eq_idx, deferred_logger); - } - - if (residual > residuals[eq_idx]) { - residuals[eq_idx] = residual; - } - } - } - - // handling the control equation residual - { - const double control_residual = std::abs(resWell_[0][numWellEq - 1]); - if (std::isnan(control_residual) || std::isinf(control_residual)) { - OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for control residal get for well " << name(), deferred_logger); - } - residuals[numWellEq] = control_residual; - } - - return residuals; - } - - - - - - template - double - MultisegmentWell:: - getResidualMeasureValue(const WellState& well_state, - const std::vector& residuals, - DeferredLogger& deferred_logger) const - { - assert(int(residuals.size()) == numWellEq + 1); - - const double rate_tolerance = param_.tolerance_wells_; - int count = 0; - double sum = 0; - for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) { - if (residuals[eq_idx] > rate_tolerance) { - sum += residuals[eq_idx] / rate_tolerance; - ++count; - } - } - - const double pressure_tolerance = param_.tolerance_pressure_ms_wells_; - if (residuals[SPres] > pressure_tolerance) { - sum += residuals[SPres] / pressure_tolerance; - ++count; - } - - const double control_tolerance = getControlTolerance(well_state, deferred_logger); - if (residuals[SPres + 1] > control_tolerance) { - sum += residuals[SPres + 1] / control_tolerance; - ++count; - } - - // if (count == 0), it should be converged. - assert(count != 0); - - return sum; - } - - - - - - template - double - MultisegmentWell:: - getControlTolerance(const WellState& well_state, - DeferredLogger& deferred_logger) const - { - double control_tolerance = 0.; - - const int well_index = index_of_well_; - if (this->isInjector() ) - { - auto current = well_state.currentInjectionControl(well_index); - switch(current) { - case Well::InjectorCMode::THP: - control_tolerance = param_.tolerance_pressure_ms_wells_; - break; - case Well::InjectorCMode::BHP: - control_tolerance = param_.tolerance_wells_; - break; - case Well::InjectorCMode::RATE: - case Well::InjectorCMode::RESV: - control_tolerance = param_.tolerance_wells_; - break; - case Well::InjectorCMode::GRUP: - control_tolerance = param_.tolerance_wells_; - break; - default: - OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); - } - } - - if (this->isProducer() ) - { - auto current = well_state.currentProductionControl(well_index); - switch(current) { - case Well::ProducerCMode::THP: - control_tolerance = param_.tolerance_pressure_ms_wells_; // 0.1 bar - break; - case Well::ProducerCMode::BHP: - control_tolerance = param_.tolerance_wells_; // 0.01 bar - break; - case Well::ProducerCMode::ORAT: - case Well::ProducerCMode::WRAT: - case Well::ProducerCMode::GRAT: - case Well::ProducerCMode::LRAT: - case Well::ProducerCMode::RESV: - case Well::ProducerCMode::CRAT: - control_tolerance = param_.tolerance_wells_; // smaller tolerance for rate control - break; - case Well::ProducerCMode::GRUP: - control_tolerance = param_.tolerance_wells_; // smaller tolerance for rate control - break; - default: - OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); - } - } - - return control_tolerance; - } - - - - - - template - void - MultisegmentWell:: - checkConvergenceControlEq(const WellState& well_state, - ConvergenceReport& report, - DeferredLogger& deferred_logger) const - { - double control_tolerance = 0.; - using CR = ConvergenceReport; - CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid; - - const int well_index = index_of_well_; - if (this->isInjector() ) - { - auto current = well_state.currentInjectionControl(well_index); - switch(current) { - case Well::InjectorCMode::THP: - ctrltype = CR::WellFailure::Type::ControlTHP; - control_tolerance = param_.tolerance_pressure_ms_wells_; - break; - case Well::InjectorCMode::BHP: - ctrltype = CR::WellFailure::Type::ControlBHP; - control_tolerance = param_.tolerance_pressure_ms_wells_; - break; - case Well::InjectorCMode::RATE: - case Well::InjectorCMode::RESV: - ctrltype = CR::WellFailure::Type::ControlRate; - control_tolerance = param_.tolerance_wells_; - break; - case Well::InjectorCMode::GRUP: - ctrltype = CR::WellFailure::Type::ControlRate; - control_tolerance = param_.tolerance_wells_; - break; - default: - OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); - } - } - - if (this->isProducer() ) - { - auto current = well_state.currentProductionControl(well_index); - switch(current) { - case Well::ProducerCMode::THP: - ctrltype = CR::WellFailure::Type::ControlTHP; - control_tolerance = param_.tolerance_pressure_ms_wells_; - break; - case Well::ProducerCMode::BHP: - ctrltype = CR::WellFailure::Type::ControlBHP; - control_tolerance = param_.tolerance_pressure_ms_wells_; - break; - case Well::ProducerCMode::ORAT: - case Well::ProducerCMode::WRAT: - case Well::ProducerCMode::GRAT: - case Well::ProducerCMode::LRAT: - case Well::ProducerCMode::RESV: - case Well::ProducerCMode::CRAT: - ctrltype = CR::WellFailure::Type::ControlRate; - control_tolerance = param_.tolerance_wells_; - break; - case Well::ProducerCMode::GRUP: - ctrltype = CR::WellFailure::Type::ControlRate; - control_tolerance = param_.tolerance_wells_; - break; - default: - OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); - } - } - - const double well_control_residual = std::abs(resWell_[0][SPres]); - const int dummy_component = -1; - const double max_residual_allowed = param_.max_residual_allowed_; - if (std::isnan(well_control_residual)) { - report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()}); - } else if (well_control_residual > max_residual_allowed * 10.) { - report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()}); - } else if ( well_control_residual > control_tolerance) { - report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()}); - } - } - - - - - - - template - void - MultisegmentWell:: - updateUpwindingSegments() - { - for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - // special treatment is needed for segment 0 - if (seg == 0) { - // we are not supposed to have injecting producers and producing injectors - assert( ! (this->isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) ); - assert( ! (this->isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) ); - upwinding_segments_[seg] = seg; - continue; - } - - // for other normal segments - if (primary_variables_evaluation_[seg][GTotal] <= 0.) { - upwinding_segments_[seg] = seg; - } else { - const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); - upwinding_segments_[seg] = outlet_segment_index; - } - } - } - - - - - template - void - MultisegmentWell:: - assembleICDPressureEq(const int seg, const UnitSystem& unit_system, - WellState& well_state, DeferredLogger& deferred_logger) const - { - // TODO: upwinding needs to be taken care of - // top segment can not be a spiral ICD device - assert(seg != 0); - - // the pressure equation is something like - // p_seg - deltaP - p_outlet = 0. - // the major part is how to calculate the deltaP - - EvalWell pressure_equation = getSegmentPressure(seg); - - EvalWell icd_pressure_drop; - switch(this->segmentSet()[seg].segmentType()) { - case Segment::SegmentType::SICD : - icd_pressure_drop = pressureDropSpiralICD(seg); - break; - case Segment::SegmentType::AICD : - icd_pressure_drop = pressureDropAutoICD(seg, unit_system); - break; - case Segment::SegmentType::VALVE : - icd_pressure_drop = pressureDropValve(seg); - break; - default: { - OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(this->segmentSet()[seg].segmentNumber()) - + " for well " + name() + " is not of ICD type", deferred_logger); - } - } - pressure_equation = pressure_equation - icd_pressure_drop; - well_state.segments(this->index_of_well_).pressure_drop_friction[seg] = icd_pressure_drop.value(); - - const int seg_upwind = upwinding_segments_[seg]; - resWell_[seg][SPres] = pressure_equation.value(); - duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq); - duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq); - } - - // contribution from the outlet segment - const int outlet_segment_index = this->segmentNumberToIndex(this->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); - - } - } - - - - template std::optional MultisegmentWell:: @@ -3109,212 +1495,6 @@ namespace Opm - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - pressureDropSpiralICD(const int seg) const - { - const SICD& sicd = this->segmentSet()[seg].spiralICD(); - - const int seg_upwind = upwinding_segments_[seg]; - const std::vector& phase_fractions = segment_phase_fractions_[seg_upwind]; - const std::vector& phase_viscosities = segment_phase_viscosities_[seg_upwind]; - - EvalWell water_fraction = 0.; - EvalWell water_viscosity = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - water_fraction = phase_fractions[water_pos]; - water_viscosity = phase_viscosities[water_pos]; - } - - EvalWell oil_fraction = 0.; - EvalWell oil_viscosity = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - oil_fraction = phase_fractions[oil_pos]; - oil_viscosity = phase_viscosities[oil_pos]; - } - - EvalWell gas_fraction = 0.; - EvalWell gas_viscosity = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - gas_fraction = phase_fractions[gas_pos]; - gas_viscosity = phase_viscosities[gas_pos]; - } - - EvalWell density = segment_densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (seg != seg_upwind) { - water_fraction.clearDerivatives(); - water_viscosity.clearDerivatives(); - oil_fraction.clearDerivatives(); - oil_viscosity.clearDerivatives(); - gas_fraction.clearDerivatives(); - gas_viscosity.clearDerivatives(); - density.clearDerivatives(); - } - - const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, - oil_fraction, oil_viscosity, sicd); - const EvalWell mixture_viscosity = (water_fraction + oil_fraction) * liquid_emulsion_viscosity + gas_fraction * gas_viscosity; - - const EvalWell reservoir_rate = segment_mass_rates_[seg] / density; - - const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor(); - - const double viscosity_cali = sicd.viscosityCalibration(); - - using MathTool = MathToolbox; - - const double density_cali = sicd.densityCalibration(); - const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75); - const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25); - - // formulation before 2016, base_strength is used - // const double base_strength = sicd.strength() / density_cali; - // formulation since 2016, strength is used instead - const double strength = sicd.strength(); - - const double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0; - - return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd; - } - - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - pressureDropAutoICD(const int seg, const UnitSystem& unit_system) const - { - const AutoICD& aicd = this->segmentSet()[seg].autoICD(); - - const int seg_upwind = this->upwinding_segments_[seg]; - const std::vector& phase_fractions = this->segment_phase_fractions_[seg_upwind]; - const std::vector& phase_viscosities = this->segment_phase_viscosities_[seg_upwind]; - const std::vector& phase_densities = this->segment_phase_densities_[seg_upwind]; - - EvalWell water_fraction = 0.; - EvalWell water_viscosity = 0.; - EvalWell water_density = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - water_fraction = phase_fractions[water_pos]; - water_viscosity = phase_viscosities[water_pos]; - water_density = phase_densities[water_pos]; - } - - EvalWell oil_fraction = 0.; - EvalWell oil_viscosity = 0.; - EvalWell oil_density = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - oil_fraction = phase_fractions[oil_pos]; - oil_viscosity = phase_viscosities[oil_pos]; - oil_density = phase_densities[oil_pos]; - } - - EvalWell gas_fraction = 0.; - EvalWell gas_viscosity = 0.; - EvalWell gas_density = 0.; - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - gas_fraction = phase_fractions[gas_pos]; - gas_viscosity = phase_viscosities[gas_pos]; - gas_density = phase_densities[gas_pos]; - } - - EvalWell density = segment_densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (seg != seg_upwind) { - water_fraction.clearDerivatives(); - water_viscosity.clearDerivatives(); - water_density.clearDerivatives(); - oil_fraction.clearDerivatives(); - oil_viscosity.clearDerivatives(); - oil_density.clearDerivatives(); - gas_fraction.clearDerivatives(); - gas_viscosity.clearDerivatives(); - gas_density.clearDerivatives(); - density.clearDerivatives(); - } - - using MathTool = MathToolbox; - const EvalWell mixture_viscosity = MathTool::pow(water_fraction, aicd.waterViscExponent()) * water_viscosity - + MathTool::pow(oil_fraction, aicd.oilViscExponent()) * oil_viscosity - + MathTool::pow(gas_fraction, aicd.gasViscExponent()) * gas_viscosity; - - const EvalWell mixture_density = MathTool::pow(water_fraction, aicd.waterDensityExponent()) * water_density - + MathTool::pow(oil_fraction, aicd.oilDensityExponent()) * oil_density - + MathTool::pow(gas_fraction, aicd.gasDensityExponent()) * gas_density; - - const double rho_reference = aicd.densityCalibration(); - const double visc_reference = aicd.viscosityCalibration(); - const auto volume_rate_icd = this->segment_mass_rates_[seg] * aicd.scalingFactor() / mixture_density; - const double sign = volume_rate_icd <= 0. ? 1.0 : -1.0; - // convert 1 unit volume rate - using M = UnitSystem::measure; - const double unit_volume_rate = unit_system.to_si(M::geometric_volume_rate, 1.); - - // TODO: we did not consider the maximum allowed rate here - const auto result = sign / rho_reference * mixture_density * mixture_density - * MathTool::pow(visc_reference/mixture_viscosity, aicd.viscExponent()) - * aicd.strength() * MathTool::pow( -sign * volume_rate_icd, aicd.flowRateExponent()) - * std::pow(unit_volume_rate, (2. - aicd.flowRateExponent())) ; - return result; - } - - - - - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - pressureDropValve(const int seg) const - { - const Valve& valve = this->segmentSet()[seg].valve(); - - const EvalWell& mass_rate = segment_mass_rates_[seg]; - const int seg_upwind = upwinding_segments_[seg]; - EvalWell visc = segment_viscosities_[seg_upwind]; - EvalWell density = segment_densities_[seg_upwind]; - // WARNING - // We disregard the derivatives from the upwind density to make sure derivatives - // wrt. to different segments dont get mixed. - if (seg != seg_upwind) { - visc.clearDerivatives(); - density.clearDerivatives(); - } - - const double additional_length = valve.pipeAdditionalLength(); - const double roughness = valve.pipeRoughness(); - const double diameter = valve.pipeDiameter(); - const double area = valve.pipeCrossArea(); - - const EvalWell friction_pressure_loss = - mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc); - - const double area_con = valve.conCrossArea(); - const double cv = valve.conFlowCoefficient(); - - const EvalWell constriction_pressure_loss = - mswellhelpers::valveContrictionPressureLoss(mass_rate, density, area_con, cv); - - const double sign = mass_rate <= 0. ? 1.0 : -1.0; - return sign * (friction_pressure_loss + constriction_pressure_loss); - } - - - - template std::vector MultisegmentWell:: @@ -3327,7 +1507,7 @@ namespace Opm const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { // calculating the perforation rate for each perforation that belongs to this segment - const EvalWell seg_pressure = getSegmentPressure(seg); + const EvalWell seg_pressure = this->getSegmentPressure(seg); for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); diff --git a/opm/simulators/wells/WellInterfaceEval.hpp b/opm/simulators/wells/WellInterfaceEval.hpp index 1a558a1d2..2557392d8 100644 --- a/opm/simulators/wells/WellInterfaceEval.hpp +++ b/opm/simulators/wells/WellInterfaceEval.hpp @@ -56,6 +56,30 @@ public: const SummaryState& summaryState, const double rho, DeferredLogger& deferred_logger) const; + template + void getGroupInjectionControl(const Group& group, + const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const InjectorType& injectorType, + const EvalWell& bhp, + const EvalWell& injection_rate, + EvalWell& control_eq, + double efficiencyFactor, + DeferredLogger& deferred_logger) const; + + + template + void getGroupProductionControl(const Group& group, + const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const EvalWell& bhp, + const std::vector& rates, + EvalWell& control_eq, + double efficiencyFactor) const; template void assembleControlEqProd(const WellState& well_state, @@ -134,55 +158,6 @@ public: protected: WellInterfaceEval(const WellInterfaceFluidSystem& baseif); - template - void getGroupInjectionControl(const Group& group, - const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const InjectorType& injectorType, - const EvalWell& bhp, - const EvalWell& injection_rate, - EvalWell& control_eq, - double efficiencyFactor, - DeferredLogger& deferred_logger) const; - - - template - void getGroupProductionControl(const Group& group, - const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const EvalWell& bhp, - const std::vector& rates, - EvalWell& control_eq, - double efficiencyFactor) const; - - template - void assembleControlEqProd(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::ProductionControls& controls, - const EvalWell& bhp, - const std::vector& rates, // Always 3 canonical rates. - const EvalWell& bhp_from_thp, - EvalWell& control_eq, - DeferredLogger& deferred_logger) const; - - template - void assembleControlEqInj_(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::InjectionControls& controls, - const EvalWell& bhp, - const EvalWell& injection_rate, - const std::function& bhp_from_thp, - EvalWell& control_eq, - DeferredLogger& deferred_logger); - const WellInterfaceFluidSystem& baseif_; }; diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index b7aae6123..b7b1dcc55 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -118,7 +118,7 @@ public: return guide_rate_; } - int numComponents() const{ + int numComponents() const { return num_components_; } diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index 284321878..7c1fc558f 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -135,11 +135,13 @@ INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,tr INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u>) // Alternative indices INSTANCE(EclAlternativeBlackOilIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>)