From eb06c4bd70209188bcad7d0a0c7ee543d8428d01 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 1 Jun 2021 15:49:24 +0200 Subject: [PATCH] add StandardWellEval --- CMakeLists_files.cmake | 1 + .../wells/MultisegmentWell_impl.hpp | 4 + opm/simulators/wells/StandardWell.hpp | 90 +- opm/simulators/wells/StandardWellEval.cpp | 1281 ++++++++++++++++ opm/simulators/wells/StandardWellEval.hpp | 195 +++ opm/simulators/wells/StandardWellGeneric.cpp | 3 +- opm/simulators/wells/StandardWellGeneric.hpp | 1 - opm/simulators/wells/StandardWell_impl.hpp | 1348 ++--------------- opm/simulators/wells/WellInterfaceEval.cpp | 12 +- opm/simulators/wells/WellInterfaceEval.hpp | 101 +- opm/simulators/wells/WellInterfaceGeneric.hpp | 8 +- opm/simulators/wells/WellInterfaceIndices.cpp | 2 + opm/simulators/wells/WellInterfaceIndices.hpp | 9 +- 13 files changed, 1689 insertions(+), 1366 deletions(-) create mode 100644 opm/simulators/wells/StandardWellEval.cpp create mode 100644 opm/simulators/wells/StandardWellEval.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index b54f4087c..74d7a150c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -66,6 +66,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/MultisegmentWellGeneric.cpp opm/simulators/wells/ParallelWellInfo.cpp opm/simulators/wells/SegmentState.cpp + opm/simulators/wells/StandardWellEval.cpp opm/simulators/wells/StandardWellGeneric.cpp opm/simulators/wells/TargetCalculator.cpp opm/simulators/wells/VFPHelpers.cpp diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 4bb07044e..2b333e631 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -27,6 +27,10 @@ #include #include +#if HAVE_CUDA || HAVE_OPENCL +#include +#endif + namespace Opm { diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index c0888dc7d..f1378b8b6 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -23,10 +23,6 @@ #ifndef OPM_STANDARDWELL_HEADER_INCLUDED #define OPM_STANDARDWELL_HEADER_INCLUDED -#if HAVE_CUDA || HAVE_OPENCL -#include -#endif - #include #include #include @@ -47,6 +43,8 @@ #include #include +#include + #include #include @@ -59,11 +57,16 @@ namespace Opm template class StandardWell : public WellInterface - , public StandardWellGeneric> + , public StandardWellEval, + GetPropType, + GetPropType> { public: typedef WellInterface Base; + using StdWellEval = StandardWellEval, + GetPropType, + GetPropType>; // TODO: some functions working with AD variables handles only with values (double) without // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value. @@ -89,6 +92,7 @@ namespace Opm using Base::has_solvent; using Base::has_zFraction; using Base::has_polymer; + using Base::has_polymermw; using Base::has_foam; using Base::has_brine; using Base::has_energy; @@ -135,9 +139,10 @@ namespace Opm using Base::Gas; using typename Base::BVector; - using typename Base::Eval; - typedef DenseAd::DynamicEvaluation EvalWell; + using Eval = typename StdWellEval::Eval; + using EvalWell = typename StdWellEval::EvalWell; + using BVectorWell = typename StdWellEval::BVectorWell; using Base::contiSolventEqIdx; using Base::contiZfracEqIdx; @@ -166,10 +171,6 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const override; - void updateWellStateWithTarget(const Simulator& ebos_simulator, - WellState& well_state, - 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, @@ -181,11 +182,6 @@ namespace Opm /// r = r - C D^-1 Rw virtual void apply(BVector& r) const override; -#if HAVE_CUDA || HAVE_OPENCL - /// add the contribution (C, D^-1, 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, @@ -329,39 +325,6 @@ namespace Opm using Base::ipr_a_; using Base::ipr_b_; using Base::changed_to_stopped_this_step_; - using typename StandardWellGeneric::BVectorWell; - - - // total number of the well equations and primary variables - // there might be extra equations be used, numWellEq will be updated during the initialization - int numWellEq_ = numStaticWellEq; - - // the values for the primary varibles - // based on different solutioin strategies, the wells can have different primary variables - mutable std::vector primary_variables_; - - // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation - mutable std::vector primary_variables_evaluation_; - - // the saturations in the well bore under surface conditions at the beginning of the time step - std::vector F0_; - - // Optimize only wells under THP control - bool glift_optimize_only_thp_wells = true; - - const EvalWell& getBhp() const; - - EvalWell getQs(const int comp_idx) const; - - const EvalWell& getWQTotal() const; - - EvalWell wellVolumeFractionScaled(const int phase) const; - - EvalWell wellVolumeFraction(const unsigned compIdx) const; - - EvalWell wellSurfaceVolumeFraction(const int phase) const; - - EvalWell extendEval(const Eval& in) const; Eval getPerfCellPressure(const FluidState& fs) const; @@ -382,14 +345,6 @@ namespace Opm std::vector& rvmax_perf, std::vector& surf_dens_perf) const; - // TODO: not total sure whether it is a good idea to put this function here - // the major reason to put here is to avoid the usage of Wells struct - void computeConnectionDensities(const std::vector& perfComponentRates, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf); - void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator, const WellState& well_state, const std::vector& b_perf, @@ -397,9 +352,6 @@ namespace Opm const std::vector& rvmax_perf, const std::vector& surf_dens_perf); - // computing the accumulation term for later use in well mass equations - void computeAccumWell(); - void computeWellConnectionPressures(const Simulator& ebosSimulator, const WellState& well_state); @@ -447,19 +399,6 @@ namespace Opm void updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const; - void updateThp(WellState& well_state, DeferredLogger& deferred_logger) const; - - - void assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - DeferredLogger& deferred_logger); - - // handle the non reasonable fractions due to numerical overshoot - void processFractions() const; - - virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator, const double dt, const Well::InjectionControls& inj_controls, @@ -508,11 +447,6 @@ namespace Opm // TODO: looking for better alternative to avoid wrong-signed well rates bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const; - // calculate a relaxation factor to avoid overshoot of the fractions for producers - // which might result in negative rates - static double relaxationFactorFractionsProducer(const std::vector& primary_variables, - const BVectorWell& dwells); - // calculate the skin pressure based on water velocity, throughput and polymer concentration. // throughput is used to describe the formation damage during water/polymer injection. // calculated skin pressure will be applied to the drawdown during perforation rate calculation diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp new file mode 100644 index 000000000..d6aaf99e7 --- /dev/null +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -0,0 +1,1281 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + Copyright 2016 - 2017 IRIS AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#if HAVE_CUDA || HAVE_OPENCL +#include +#endif + + +namespace Opm +{ + +template +StandardWellEval:: +StandardWellEval(const WellInterfaceIndices& baseif) + : StandardWellGeneric(Bhp, baseif) + , baseif_(baseif) + , F0_(numWellConservationEq) +{ +} + +template +void StandardWellEval:: +initPrimaryVariablesEvaluation() const +{ + for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) { + primary_variables_evaluation_[eqIdx] = + EvalWell::createVariable(numWellEq_ + Indices::numEq, primary_variables_[eqIdx], Indices::numEq + eqIdx); + } +} + +template +typename StandardWellEval::EvalWell +StandardWellEval:: +extendEval(const Eval& in) const +{ + EvalWell out(numWellEq_ + Indices::numEq, in.value()); + for(int eqIdx = 0; eqIdx < Indices::numEq;++eqIdx) { + out.setDerivative(eqIdx, in.derivative(eqIdx)); + } + return out; +} + +template +double +StandardWellEval:: +relaxationFactorFractionsProducer(const std::vector& primary_variables, + const BVectorWell& dwells) +{ + // TODO: not considering solvent yet + // 0.95 is a experimental value, which remains to be optimized + double relaxation_factor = 1.0; + + if (FluidSystem::numActivePhases() > 1) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const double relaxation_factor_w = StandardWellGeneric:: + relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const double relaxation_factor_g = StandardWellGeneric:: + relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will + // not be negative oil fraction later + const double original_sum = primary_variables[WFrac] + primary_variables[GFrac]; + const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor; + const double possible_updated_sum = original_sum - relaxed_update; + + if (possible_updated_sum > 1.0) { + assert(relaxed_update != 0.); + + const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95; + relaxation_factor *= further_relaxation_factor; + } + } + + assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0); + } + return relaxation_factor; +} + +template +typename StandardWellEval::EvalWell +StandardWellEval:: +wellVolumeFraction(const unsigned compIdx) const +{ + if (FluidSystem::numActivePhases() == 1) { + return EvalWell(numWellEq_ + Indices::numEq, 1.0); + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { + return primary_variables_evaluation_[WFrac]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { + return primary_variables_evaluation_[GFrac]; + } + + if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) { + return primary_variables_evaluation_[SFrac]; + } + } + else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { + return primary_variables_evaluation_[GFrac]; + } + } + + // Oil or WATER fraction + EvalWell well_fraction(numWellEq_ + Indices::numEq, 1.0); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + well_fraction -= primary_variables_evaluation_[WFrac]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + well_fraction -= primary_variables_evaluation_[GFrac]; + } + + if (Indices::enableSolvent) { + well_fraction -= primary_variables_evaluation_[SFrac]; + } + } + else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))) { + + well_fraction -= primary_variables_evaluation_[GFrac]; + } + + return well_fraction; +} + +template +typename StandardWellEval::EvalWell +StandardWellEval:: +getQs(const int comp_idx) const +{ + // Note: currently, the WQTotal definition is still depends on Injector/Producer. + assert(comp_idx < baseif_.numComponents()); + + if (baseif_.isInjector()) { // only single phase injection + double inj_frac = 0.0; + switch (baseif_.wellEcl().injectorType()) { + case InjectorType::WATER: + if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) { + inj_frac = 1.0; + } + break; + case InjectorType::GAS: + if (Indices::enableSolvent && comp_idx == Indices::contiSolventEqIdx) { // solvent + inj_frac = baseif_.wsolvent(); + } else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) { + inj_frac = Indices::enableSolvent ? 1.0 - baseif_.wsolvent() : 1.0; + } + break; + case InjectorType::OIL: + if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) { + inj_frac = 1.0; + } + break; + case InjectorType::MULTI: + // Not supported. + // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + // "Multi phase injectors are not supported, requested for well " + name()); + break; + } + return inj_frac * primary_variables_evaluation_[WQTotal]; + } else { // producers + return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx); + } +} + +template +typename StandardWellEval::EvalWell +StandardWellEval:: +wellVolumeFractionScaled(const int compIdx) const +{ + const int legacyCompIdx = baseif_.ebosCompIdxToFlowCompIdx(compIdx); + const double scal = baseif_.scalingFactor(legacyCompIdx); + if (scal > 0) + return wellVolumeFraction(compIdx) / scal; + + // the scaling factor may be zero for RESV controlled wells. + return wellVolumeFraction(compIdx); +} + +template +typename StandardWellEval::EvalWell +StandardWellEval:: +wellSurfaceVolumeFraction(const int compIdx) const +{ + EvalWell sum_volume_fraction_scaled(numWellEq_ + Indices::numEq, 0.); + for (int idx = 0; idx < baseif_.numComponents(); ++idx) { + sum_volume_fraction_scaled += wellVolumeFractionScaled(idx); + } + + assert(sum_volume_fraction_scaled.value() != 0.); + + return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled; + } + +template +void +StandardWellEval:: +updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = WellInterfaceIndices::Gas; + static constexpr int Oil = WellInterfaceIndices::Oil; + static constexpr int Water = WellInterfaceIndices::Water; + + if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return; + + const int well_index = baseif_.indexOfWell(); + const int np = baseif_.numPhases(); + const auto& pu = baseif_.phaseUsage(); + + // the weighted total well rate + double total_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + total_well_rate += baseif_.scalingFactor(p) * well_state.wellRates(well_index)[p]; + } + + // Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate + // under surface condition is used here + if (baseif_.isInjector()) { + switch (baseif_.wellEcl().injectorType()) { + case InjectorType::WATER: + primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Water]]; + break; + case InjectorType::GAS: + primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Gas]]; + break; + case InjectorType::OIL: + primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Oil]]; + break; + case InjectorType::MULTI: + // Not supported. + deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + "Multi phase injectors are not supported, requested for well " + baseif_.name()); + break; + } + } else { + primary_variables_[WQTotal] = total_well_rate; + } + + if (std::abs(total_well_rate) > 0.) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * well_state.wellRates(well_index)[pu.phase_pos[Water]] / total_well_rate; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates(well_index)[pu.phase_pos[Gas]] + - (Indices::enableSolvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ; + } + if constexpr (Indices::enableSolvent) { + primary_variables_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; + } + } else { // total_well_rate == 0 + if (baseif_.isInjector()) { + auto phase = baseif_.wellEcl().getInjectionProperties().injectorType; + // only single phase injection handled + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (phase == InjectorType::WATER) { + primary_variables_[WFrac] = 1.0; + } else { + primary_variables_[WFrac] = 0.0; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + if (phase == InjectorType::GAS) { + primary_variables_[GFrac] = 1.0; + if constexpr (Indices::enableSolvent) { + primary_variables_[GFrac] = 1.0 - baseif_.wsolvent(); + primary_variables_[SFrac] = baseif_.wsolvent(); + } + } else { + primary_variables_[GFrac] = 0.0; + } + } + + // TODO: it is possible to leave injector as a oil well, + // when F_w and F_g both equals to zero, not sure under what kind of circumstance + // this will happen. + } else if (baseif_.isProducer()) { // producers + // TODO: the following are not addressed for the solvent case yet + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + primary_variables_[WFrac] = 1.0 / np; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + primary_variables_[GFrac] = 1.0 / np; + } + } else { + OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger); + } + } + + + // BHP + primary_variables_[Bhp] = well_state.bhp(baseif_.indexOfWell()); +} + +template +void +StandardWellEval:: +assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + DeferredLogger& deferred_logger) +{ + static constexpr int Gas = WellInterfaceIndices::Gas; + static constexpr int Oil = WellInterfaceIndices::Oil; + static constexpr int Water = WellInterfaceIndices::Water; + EvalWell control_eq(numWellEq_ + Indices::numEq, 0.0); + + const auto& well = baseif_.wellEcl(); + + auto getRates = [&]() { + std::vector rates(3, EvalWell(numWellEq_ + Indices::numEq, 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 injection rate. + const EvalWell injection_rate = getWQTotal(); + // 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, this->getRho(), deferred_logger); + }; + // Call generic implementation. + const auto& inj_controls = well.injectionControls(summaryState); + 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, this->getRho(), deferred_logger); + }; + // Call generic implementation. + const auto& prod_controls = well.productionControls(summaryState); + 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 + // TODO: we should use a different index system for the well equations + this->resWell_[0][Bhp] = control_eq.value(); + for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) { + this->invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); + } +} + + +template +void +StandardWellEval:: +updatePrimaryVariablesPolyMW(const BVectorWell& dwells) const +{ + if (baseif_.isInjector()) { + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + const int wat_vel_index = Bhp + 1 + perf; + const int pskin_index = Bhp + 1 + baseif_.numPerfs() + perf; + + const double relaxation_factor = 0.9; + const double dx_wat_vel = dwells[0][wat_vel_index]; + primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel; + + const double dx_pskin = dwells[0][pskin_index]; + primary_variables_[pskin_index] -= relaxation_factor * dx_pskin; + } + } +} + +template +void +StandardWellEval:: +processFractions() const +{ + static constexpr int Gas = WellInterfaceIndices::Gas; + static constexpr int Oil = WellInterfaceIndices::Oil; + static constexpr int Water = WellInterfaceIndices::Water; + const auto pu = baseif_.phaseUsage(); + std::vector F(baseif_.numPhases(), 0.0); + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + F[pu.phase_pos[Oil]] = 1.0; + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + F[pu.phase_pos[Water]] = primary_variables_[WFrac]; + F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; + F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]]; + } + } + else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + F[pu.phase_pos[Water]] = 1.0; + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; + F[pu.phase_pos[Water]] -= F[pu.phase_pos[Gas]]; + } + } + else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + F[pu.phase_pos[Gas]] = 1.0; + } + + [[maybe_unused]] double F_solvent; + if constexpr (Indices::enableSolvent) { + F_solvent = primary_variables_[SFrac]; + F[pu.phase_pos[Oil]] -= F_solvent; + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (F[Water] < 0.0) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]); + } + if constexpr (Indices::enableSolvent) { + F_solvent /= (1.0 - F[pu.phase_pos[Water]]); + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]); + } + F[pu.phase_pos[Water]] = 0.0; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + if (F[pu.phase_pos[Gas]] < 0.0) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); + } + if constexpr (Indices::enableSolvent) { + F_solvent /= (1.0 - F[pu.phase_pos[Gas]]); + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]); + } + F[pu.phase_pos[Gas]] = 0.0; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + if (F[pu.phase_pos[Oil]] < 0.0) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); + } + if constexpr (Indices::enableSolvent) { + F_solvent /= (1.0 - F[pu.phase_pos[Oil]]); + } + F[pu.phase_pos[Oil]] = 0.0; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + primary_variables_[WFrac] = F[pu.phase_pos[Water]]; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; + } + if constexpr (Indices::enableSolvent) { + primary_variables_[SFrac] = F_solvent; + } +} + + +template +void +StandardWellEval:: +updateThp(WellState& well_state, + DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = WellInterfaceIndices::Gas; + static constexpr int Oil = WellInterfaceIndices::Oil; + static constexpr int Water = WellInterfaceIndices::Water; + + // 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(well_state, + rates, + bhp, + deferred_logger)); +} + +template +void +StandardWellEval:: +updateWellStateFromPrimaryVariables(WellState& well_state, + DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = WellInterfaceIndices::Gas; + static constexpr int Oil = WellInterfaceIndices::Oil; + static constexpr int Water = WellInterfaceIndices::Water; + + const PhaseUsage& pu = baseif_.phaseUsage(); + std::vector F(baseif_.numPhases(), 0.0); + [[maybe_unused]] double F_solvent = 0.0; + if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) { + const int oil_pos = pu.phase_pos[Oil]; + F[oil_pos] = 1.0; + + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int water_pos = pu.phase_pos[Water]; + F[water_pos] = primary_variables_[WFrac]; + F[oil_pos] -= F[water_pos]; + } + + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int gas_pos = pu.phase_pos[Gas]; + F[gas_pos] = primary_variables_[GFrac]; + F[oil_pos] -= F[gas_pos]; + } + + if constexpr (Indices::enableSolvent) { + F_solvent = primary_variables_[SFrac]; + F[oil_pos] -= F_solvent; + } + } + else if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int water_pos = pu.phase_pos[Water]; + F[water_pos] = 1.0; + + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int gas_pos = pu.phase_pos[Gas]; + F[gas_pos] = primary_variables_[GFrac]; + F[water_pos] -= F[gas_pos]; + } + } + else if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int gas_pos = pu.phase_pos[Gas]; + F[gas_pos] = 1.0; + } + + + // convert the fractions to be Q_p / G_total to calculate the phase rates + for (int p = 0; p < baseif_.numPhases(); ++p) { + const double scal = baseif_.scalingFactor(p); + // for injection wells, there should only one non-zero scaling factor + if (scal > 0) { + F[p] /= scal ; + } else { + // this should only happens to injection wells + F[p] = 0.; + } + } + + // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. + // More testing is needed to make sure this is correct for well groups and THP. + if constexpr (Indices::enableSolvent){ + F_solvent /= baseif_.scalingFactor(Indices::contiSolventEqIdx); + F[pu.phase_pos[Gas]] += F_solvent; + } + + well_state.update_bhp(baseif_.indexOfWell(), primary_variables_[Bhp]); + + // calculate the phase rates based on the primary variables + // for producers, this is not a problem, while not sure for injectors here + if (baseif_.isProducer()) { + const double g_total = primary_variables_[WQTotal]; + for (int p = 0; p < baseif_.numPhases(); ++p) { + well_state.wellRates(baseif_.indexOfWell())[p] = g_total * F[p]; + } + } else { // injectors + for (int p = 0; p < baseif_.numPhases(); ++p) { + well_state.wellRates(baseif_.indexOfWell())[p] = 0.0; + } + switch (baseif_.wellEcl().injectorType()) { + case InjectorType::WATER: + well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Water]] = primary_variables_[WQTotal]; + break; + case InjectorType::GAS: + well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Gas]] = primary_variables_[WQTotal]; + break; + case InjectorType::OIL: + well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Oil]] = primary_variables_[WQTotal]; + break; + case InjectorType::MULTI: + // Not supported. + deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + "Multi phase injectors are not supported, requested for well " + baseif_.name()); + break; + } + } + + updateThp(well_state, deferred_logger); +} + +template +void +StandardWellEval:: +updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const +{ + if (baseif_.isInjector()) { + auto& perf_water_velocity = well_state.perfWaterVelocity(baseif_.indexOfWell()); + auto& perf_skin_pressure = well_state.perfSkinPressure(baseif_.indexOfWell()); + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + perf_water_velocity[perf] = primary_variables_[Bhp + 1 + perf]; + perf_skin_pressure[perf] = primary_variables_[Bhp + 1 + baseif_.numPerfs() + perf]; + } + } +} + +template +void +StandardWellEval:: +computeAccumWell() +{ + for (int eq_idx = 0; eq_idx < numWellConservationEq; ++eq_idx) { + F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value(); + } +} + +template +void +StandardWellEval:: +updatePrimaryVariablesNewton(const BVectorWell& dwells, + const double dFLimit, + const double dBHPLimit) const +{ + const std::vector old_primary_variables = primary_variables_; + + // for injectors, very typical one of the fractions will be one, and it is easy to get zero value + // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here + const double relaxation_factor_fractions = + (baseif_.isProducer()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells) + : 1.0; + + // update the second and third well variable (The flux fractions) + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); + // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; + primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); + primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; + } + + if constexpr (Indices::enableSolvent) { + const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); + primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; + } + + processFractions(); + + // updating the total rates Q_t + const double relaxation_factor_rate = this->relaxationFactorRate(old_primary_variables, dwells); + primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate; + + // updating the bottom hole pressure + { + const int sign1 = dwells[0][Bhp] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit); + // 1e5 to make sure bhp will not be below 1bar + primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5); + } +} + +template +ConvergenceReport +StandardWellEval:: +getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + const double tol_wells, + const double maxResidualAllowed, + std::vector& res, + DeferredLogger& deferred_logger) const +{ + res.resize(numWellEq_); + for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) { + // magnitude of the residual matters + res[eq_idx] = std::abs(this->resWell_[0][eq_idx]); + } + + std::vector well_flux_residual(baseif_.numComponents()); + + // Finish computation + for (int compIdx = 0; compIdx < baseif_.numComponents(); ++compIdx ) + { + well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx]; + } + + ConvergenceReport report; + using CR = ConvergenceReport; + CR::WellFailure::Type type = CR::WellFailure::Type::MassBalance; + // checking if any NaN or too large residuals found + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + const int compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); + + if (std::isnan(well_flux_residual[compIdx])) { + report.setWellFailed({type, CR::Severity::NotANumber, compIdx, baseif_.name()}); + } else if (well_flux_residual[compIdx] > maxResidualAllowed) { + report.setWellFailed({type, CR::Severity::TooLarge, compIdx, baseif_.name()}); + } else if (well_flux_residual[compIdx] > tol_wells) { + report.setWellFailed({type, CR::Severity::Normal, compIdx, baseif_.name()}); + } + } + + this->checkConvergenceControlEq(well_state, report, deferred_logger, maxResidualAllowed); + + return report; +} + +template +void +StandardWellEval:: +computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf) +{ + // Verify that we have consistent input. + const int nperf = baseif_.numPerfs(); + const int num_comp = baseif_.numComponents(); + + // 1. Compute the flow (in surface volume units for each + // component) exiting up the wellbore from each perforation, + // taking into account flow from lower in the well, and + // in/out-flow at each perforation. + std::vector q_out_perf((nperf)*num_comp, 0.0); + + // Step 1 depends on the order of the perforations. Hence we need to + // do the modifications globally. + // Create and get the global perforation information and do this sequentially + // on each process + + const auto& factory = baseif_.parallelWellInfo().getGlobalPerfContainerFactory(); + auto global_q_out_perf = factory.createGlobal(q_out_perf, num_comp); + auto global_perf_comp_rates = factory.createGlobal(perfComponentRates, num_comp); + + // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore + // Iterate over well perforations from bottom to top. + for (int perf = factory.numGlobalPerfs() - 1; perf >= 0; --perf) { + for (int component = 0; component < num_comp; ++component) { + auto index = perf * num_comp + component; + if (perf == factory.numGlobalPerfs() - 1) { + // This is the bottom perforation. No flow from below. + global_q_out_perf[index] = 0.0; + } else { + // Set equal to flow from below. + global_q_out_perf[index] = global_q_out_perf[index + num_comp]; + } + // Subtract outflow through perforation. + global_q_out_perf[index] -= global_perf_comp_rates[index]; + } + } + + // Copy the data back to local view + factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_comp); + + // 2. Compute the component mix at each perforation as the + // absolute values of the surface rates divided by their sum. + // Then compute volume ratios (formation factors) for each perforation. + // Finally compute densities for the segments associated with each perforation. + std::vector mix(num_comp,0.0); + std::vector x(num_comp); + std::vector surf_dens(num_comp); + + for (int perf = 0; perf < nperf; ++perf) { + // Find component mix. + const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf, + q_out_perf.begin() + num_comp*(perf+1), 0.0); + if (tot_surf_rate != 0.0) { + for (int component = 0; component < num_comp; ++component) { + mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); + } + } else if (num_comp == 1) { + mix[num_comp-1] = 1.0; + } else { + std::fill(mix.begin(), mix.end(), 0.0); + // No flow => use well specified fractions for mix. + if (baseif_.isInjector()) { + switch (baseif_.wellEcl().injectorType()) { + case InjectorType::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + case InjectorType::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case InjectorType::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case InjectorType::MULTI: + // Not supported. + // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + // "Multi phase injectors are not supported, requested for well " + name()); + break; + } + } else { + assert(baseif_.isProducer()); + // For the frist perforation without flow we use the preferred phase to decide the mix initialization. + if (perf == 0) { // + switch (baseif_.wellEcl().getPreferredPhase()) { + case Phase::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case Phase::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case Phase::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + default: + // No others supported. + break; + } + // For the rest of the perforation without flow we use mix from the above perforation. + } else { + mix = x; + } + + } + } + // Compute volume ratio. + x = mix; + + // Subtract dissolved gas from oil phase and vapporized oil from gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + double rs = 0.0; + double rv = 0.0; + if (!rsmax_perf.empty() && mix[oilpos] > 1e-12) { + rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); + } + if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { + rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); + } + if (rs != 0.0) { + // Subtract gas in oil from gas mixture + x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv); + } + if (rv != 0.0) { + // Subtract oil in gas from oil mixture + x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);; + } + } + double volrat = 0.0; + for (int component = 0; component < num_comp; ++component) { + volrat += x[component] / b_perf[perf*num_comp+ component]; + } + for (int component = 0; component < num_comp; ++component) { + surf_dens[component] = surf_dens_perf[perf*num_comp+ component]; + } + + // Compute segment density. + this->perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; + } +} + +template +void +StandardWellEval:: +computePerfRate(const std::vector& mob, + const EvalWell& pressure, + const EvalWell& bhp, + const EvalWell& rs, + const EvalWell& rv, + std::vector& b_perfcells_dense, + const double Tw, + const int perf, + const bool allow_cf, + const bool enable_polymermw, + std::vector& cq_s, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const +{ + // Pressure drawdown (also used to determine direction of flow) + const EvalWell well_pressure = bhp + this->perf_pressure_diffs_[perf]; + EvalWell drawdown = pressure - well_pressure; + + if (enable_polymermw) { + if (baseif_.isInjector()) { + const int pskin_index = Bhp + 1 + baseif_.numPerfs() + perf; + const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index]; + drawdown += skin_pressure; + } + } + + // producing perforations + if ( drawdown.value() > 0 ) { + //Do nothing if crossflow is not allowed + if (!allow_cf && baseif_.isInjector()) { + return; + } + + // compute component volumetric rates at standard conditions + for (int componentIdx = 0; componentIdx < baseif_.numComponents(); ++componentIdx) { + const EvalWell cq_p = - Tw * (mob[componentIdx] * drawdown); + cq_s[componentIdx] = b_perfcells_dense[componentIdx] * 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_sOil = cq_s[oilCompIdx]; + const EvalWell cq_sGas = cq_s[gasCompIdx]; + const EvalWell dis_gas = rs * cq_sOil; + const EvalWell vap_oil = rv * cq_sGas; + + cq_s[gasCompIdx] += dis_gas; + cq_s[oilCompIdx] += vap_oil; + + // recording the perforation solution gas rate and solution oil rates + if (baseif_.isProducer()) { + perf_dis_gas_rate = dis_gas.value(); + perf_vap_oil_rate = vap_oil.value(); + } + } + + } else { + //Do nothing if crossflow is not allowed + if (!allow_cf && baseif_.isProducer()) { + return; + } + + // Using total mobilities + EvalWell total_mob_dense = mob[0]; + for (int componentIdx = 1; componentIdx < baseif_.numComponents(); ++componentIdx) { + total_mob_dense += mob[componentIdx]; + } + + // injection perforations total volume rates + const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); + + // surface volume fraction of fluids within wellbore + std::vector cmix_s(baseif_.numComponents(), EvalWell{numWellEq_ + Indices::numEq}); + for (int componentIdx = 0; componentIdx < baseif_.numComponents(); ++componentIdx) { + cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); + } + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio(numWellEq_ + Indices::numEq, 0.); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; + } + + if constexpr (Indices::enableSolvent) { + volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx]; + } + + 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 + const EvalWell d = EvalWell(numWellEq_ + Indices::numEq, 1.0) - rv * rs; + + if (d.value() == 0.0) { + OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << baseif_.name() << " during flux calcuation" + << " with rs " << rs << " and rv " << rv, deferred_logger); + } + + const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; + //std::cout << "tmp_oil " < +void +StandardWellEval:: +init(std::vector& perf_depth, + const std::vector& depth_arg, + const int num_cells, + const bool has_polymermw) +{ + perf_depth.resize(baseif_.numPerfs(), 0.); + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + const int cell_idx = baseif_.cells()[perf]; + perf_depth[perf] = depth_arg[cell_idx]; + } + + // counting/updating primary variable numbers + if (has_polymermw) { + if (baseif_.isInjector()) { + // adding a primary variable for water perforation rate per connection + numWellEq_ += baseif_.numPerfs(); + // adding a primary variable for skin pressure per connection + numWellEq_ += baseif_.numPerfs(); + } + } + + // with the updated numWellEq_, we can initialize the primary variables and matrices now + primary_variables_.resize(numWellEq_, 0.0); + primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + Indices::numEq, 0.0}); + + // setup sparsity pattern for the matrices + //[A C^T [x = [ res + // B D] x_well] res_well] + // set the size of the matrices + this->invDuneD_.setSize(1, 1, 1); + this->duneB_.setSize(1, num_cells, baseif_.numPerfs()); + this->duneC_.setSize(1, num_cells, baseif_.numPerfs()); + + for (auto row=this->invDuneD_.createbegin(), end = this->invDuneD_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + row.insert(row.index()); + } + // the block size is run-time determined now + this->invDuneD_[0][0].resize(numWellEq_, numWellEq_); + + for (auto row = this->duneB_.createbegin(), end = this->duneB_.createend(); row!=end; ++row) { + for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) { + const int cell_idx = baseif_.cells()[perf]; + row.insert(cell_idx); + } + } + + for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) { + const int cell_idx = baseif_.cells()[perf]; + // the block size is run-time determined now + this->duneB_[0][cell_idx].resize(numWellEq_, Indices::numEq); + } + + // make the C^T matrix + for (auto row = this->duneC_.createbegin(), end = this->duneC_.createend(); row!=end; ++row) { + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + const int cell_idx = baseif_.cells()[perf]; + row.insert(cell_idx); + } + } + + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + const int cell_idx = baseif_.cells()[perf]; + this->duneC_[0][cell_idx].resize(numWellEq_, Indices::numEq); + } + + this->resWell_.resize(1); + // the block size of resWell_ is also run-time determined now + this->resWell_[0].resize(numWellEq_); + + // resize temporary class variables + this->Bx_.resize( this->duneB_.N() ); + for (unsigned i = 0; i < this->duneB_.N(); ++i) { + this->Bx_[i].resize(numWellEq_); + } + + this->invDrw_.resize( this->invDuneD_.N() ); + for (unsigned i = 0; i < this->invDuneD_.N(); ++i) { + this->invDrw_[i].resize(numWellEq_); + } +} + +#if HAVE_CUDA || HAVE_OPENCL +template +void +StandardWellEval:: +addWellContribution(WellContributions& wellContribs) const +{ + std::vector colIndices; + std::vector nnzValues; + colIndices.reserve(this->duneB_.nonzeroes()); + nnzValues.reserve(this->duneB_.nonzeroes()*numStaticWellEq * Indices::numEq); + + // duneC + for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC ) + { + colIndices.emplace_back(colC.index()); + for (int i = 0; i < numStaticWellEq; ++i) { + for (int j = 0; j < Indices::numEq; ++j) { + nnzValues.emplace_back((*colC)[i][j]); + } + } + } + wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->duneC_.nonzeroes()); + + // invDuneD + colIndices.clear(); + nnzValues.clear(); + colIndices.emplace_back(0); + for (int i = 0; i < numStaticWellEq; ++i) + { + for (int j = 0; j < numStaticWellEq; ++j) { + nnzValues.emplace_back(this->invDuneD_[0][0][i][j]); + } + } + wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1); + + // duneB + colIndices.clear(); + nnzValues.clear(); + for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB ) + { + colIndices.emplace_back(colB.index()); + for (int i = 0; i < numStaticWellEq; ++i) { + for (int j = 0; j < Indices::numEq; ++j) { + nnzValues.emplace_back((*colB)[i][j]); + } + } + } + wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->duneB_.nonzeroes()); +} +#endif + +#define INSTANCE(A,...) \ +template class StandardWellEval,__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,1u,0u,false,true,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<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>) + +// Alternative indices +INSTANCE(EclAlternativeBlackOilIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) + +} diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp new file mode 100644 index 000000000..6ba3ca3dd --- /dev/null +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -0,0 +1,195 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + Copyright 2016 - 2017 IRIS AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_STANDARDWELL_EVAL_HEADER_INCLUDED +#define OPM_STANDARDWELL_EVAL_HEADER_INCLUDED + +#include + +#include + +#include +#include + +namespace Opm +{ + +class ConvergenceReport; +class DeferredLogger; +class GroupState; +class Schedule; +class SummaryState; +class WellContributions; +template class WellInterfaceIndices; +class WellState; + +template +class StandardWellEval : public StandardWellGeneric +{ +protected: + // number of the conservation equations + static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents; + // number of the well control equations + static constexpr int numWellControlEq = 1; + // number of the well equations that will always be used + // based on the solution strategy, there might be other well equations be introduced + static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq; + // the index for Bhp in primary variables and also the index of well control equation + // they both will be the last one in their respective system. + // TODO: we should have indices for the well equations and well primary variables separately + static constexpr int Bhp = numStaticWellEq - numWellControlEq; + + // the positions of the primary variables for StandardWell + // the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g, + // which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP. + // correspondingly, we have four well equations for blackoil model, the first three are mass + // converstation equations, and the last one is the well control equation. + // primary variables related to other components, will be before the Bhp and after F_g. + // well control equation is always the last well equation. + // TODO: in the current implementation, we use the well rate as the first primary variables for injectors, + // instead of G_t. + static constexpr bool gasoil = Indices::numPhases == 2 && (Indices::compositionSwitchIdx >= 0); + static constexpr int WQTotal = 0; + static constexpr int WFrac = gasoil ? -1000 : 1; + static constexpr int GFrac = gasoil ? 1 : 2; + static constexpr int SFrac = !Indices::enableSolvent ? -1000 : 3; + +public: + using EvalWell = DenseAd::DynamicEvaluation; + using Eval = DenseAd::Evaluation; + using BVectorWell = typename StandardWellGeneric::BVectorWell; + +#if HAVE_CUDA || HAVE_OPENCL + /// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object + void addWellContribution(WellContributions& wellContribs) const; +#endif + +protected: + StandardWellEval(const WellInterfaceIndices& baseif); + + const WellInterfaceIndices& baseif_; + + void initPrimaryVariablesEvaluation() const; + + const EvalWell& getBhp() const + { + return primary_variables_evaluation_[Bhp]; + } + + const EvalWell& getWQTotal() const + { + return primary_variables_evaluation_[WQTotal]; + } + + EvalWell extendEval(const Eval& in) const; + EvalWell getQs(const int compIdx) const; + EvalWell wellSurfaceVolumeFraction(const int compIdx) const; + EvalWell wellVolumeFraction(const unsigned compIdx) const; + EvalWell wellVolumeFractionScaled(const int phase) const; + + // calculate a relaxation factor to avoid overshoot of the fractions for producers + // which might result in negative rates + static double relaxationFactorFractionsProducer(const std::vector& primary_variables, + const BVectorWell& dwells); + + void assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + DeferredLogger& deferred_logger); + + // computing the accumulation term for later use in well mass equations + void computeAccumWell(); + + // TODO: not total sure whether it is a good idea to put this function here + // the major reason to put here is to avoid the usage of Wells struct + void computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf); + + void computePerfRate(const std::vector& mob, + const EvalWell& pressure, + const EvalWell& bhp, + const EvalWell& rs, + const EvalWell& rv, + std::vector& b_perfcells_dense, + const double Tw, + const int perf, + const bool allow_cf, + const bool enable_polymermw, + std::vector& cq_s, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const; + + ConvergenceReport getWellConvergence(const WellState& well_state, + const std::vector& B_avg, + const double tol_wells, + const double maxResidualAllowed, + std::vector& res, + DeferredLogger& deferred_logger) const; + + void init(std::vector& perf_depth, + const std::vector& depth_arg, + const int num_cells, + const bool has_polymermw); + + // handle the non reasonable fractions due to numerical overshoot + void processFractions() const; + + void updatePrimaryVariables(const WellState& well_state, + DeferredLogger& deferred_logger) const; + + void updatePrimaryVariablesPolyMW(const BVectorWell& dwells) const; + + void updateWellStateFromPrimaryVariables(WellState& well_state, + DeferredLogger& deferred_logger) const; + + void updatePrimaryVariablesNewton(const BVectorWell& dwells, + const double dFLimit, + const double dBHPLimit) const; + + void updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const; + + void updateThp(WellState& well_state, + DeferredLogger& deferred_logger) const; + + // total number of the well equations and primary variables + // there might be extra equations be used, numWellEq will be updated during the initialization + int numWellEq_ = numStaticWellEq; + + // the values for the primary varibles + // based on different solutioin strategies, the wells can have different primary variables + mutable std::vector primary_variables_; + + // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation + mutable std::vector primary_variables_evaluation_; + + // the saturations in the well bore under surface conditions at the beginning of the time step + std::vector F0_; +}; + +} + +#endif // OPM_STANDARDWELL_EVAL_HEADER_INCLUDED diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp index 553810d2d..022cb18cf 100644 --- a/opm/simulators/wells/StandardWellGeneric.cpp +++ b/opm/simulators/wells/StandardWellGeneric.cpp @@ -48,12 +48,11 @@ namespace Opm template StandardWellGeneric:: StandardWellGeneric(int Bhp, - const ParallelWellInfo& pw_info, const WellInterfaceGeneric& baseif) : baseif_(baseif) , perf_densities_(baseif_.numPerfs()) , perf_pressure_diffs_(baseif_.numPerfs()) - , parallelB_(duneB_, pw_info) + , parallelB_(duneB_, baseif_.parallelWellInfo()) , Bhp_(Bhp) { duneB_.setBuildMode(OffDiagMatWell::row_wise); diff --git a/opm/simulators/wells/StandardWellGeneric.hpp b/opm/simulators/wells/StandardWellGeneric.hpp index abd161b96..057cd01b9 100644 --- a/opm/simulators/wells/StandardWellGeneric.hpp +++ b/opm/simulators/wells/StandardWellGeneric.hpp @@ -72,7 +72,6 @@ public: protected: StandardWellGeneric(int Bhp, - const ParallelWellInfo& pw_info, const WellInterfaceGeneric& baseif); // calculate a relaxation factor to avoid overshoot of total rates diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index f21352c1c..e583a6ab1 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -45,8 +45,7 @@ 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) - , StandardWellGeneric(Bhp, pw_info, *this) - , F0_(numWellConservationEq) + , StdWellEval(static_cast&>(*this)) { assert(num_components_ == numWellConservationEq); } @@ -65,82 +64,7 @@ namespace Opm const std::vector< Scalar >& B_avg) { Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg); - - perf_depth_.resize(number_of_perforations_, 0.); - for (int perf = 0; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - perf_depth_[perf] = depth_arg[cell_idx]; - } - - // counting/updating primary variable numbers - if constexpr (Base::has_polymermw) { - if (this->isInjector()) { - // adding a primary variable for water perforation rate per connection - numWellEq_ += number_of_perforations_; - // adding a primary variable for skin pressure per connection - numWellEq_ += number_of_perforations_; - } - } - - // with the updated numWellEq_, we can initialize the primary variables and matrices now - primary_variables_.resize(numWellEq_, 0.0); - primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + numEq, 0.0}); - - // setup sparsity pattern for the matrices - //[A C^T [x = [ res - // B D] x_well] res_well] - // set the size of the matrices - this->invDuneD_.setSize(1, 1, 1); - this->duneB_.setSize(1, num_cells, number_of_perforations_); - this->duneC_.setSize(1, num_cells, number_of_perforations_); - - for (auto row=this->invDuneD_.createbegin(), end = this->invDuneD_.createend(); row!=end; ++row) { - // Add nonzeros for diagonal - row.insert(row.index()); - } - // the block size is run-time determined now - this->invDuneD_[0][0].resize(numWellEq_, numWellEq_); - - for (auto row = this->duneB_.createbegin(), end = this->duneB_.createend(); row!=end; ++row) { - for (int perf = 0 ; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - row.insert(cell_idx); - } - } - - for (int perf = 0 ; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - // the block size is run-time determined now - this->duneB_[0][cell_idx].resize(numWellEq_, numEq); - } - - // make the C^T matrix - for (auto row = this->duneC_.createbegin(), end = this->duneC_.createend(); row!=end; ++row) { - for (int perf = 0; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - row.insert(cell_idx); - } - } - - for (int perf = 0; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - this->duneC_[0][cell_idx].resize(numWellEq_, numEq); - } - - this->resWell_.resize(1); - // the block size of resWell_ is also run-time determined now - this->resWell_[0].resize(numWellEq_); - - // resize temporary class variables - this->Bx_.resize( this->duneB_.N() ); - for (unsigned i = 0; i < this->duneB_.N(); ++i) { - this->Bx_[i].resize(numWellEq_); - } - - this->invDrw_.resize( this->invDuneD_.N() ); - for (unsigned i = 0; i < this->invDuneD_.N(); ++i) { - this->invDrw_[i].resize(numWellEq_); - } + this->StdWellEval::init(perf_depth_, depth_arg, num_cells, Base::has_polymermw); } @@ -151,189 +75,7 @@ namespace Opm void StandardWell:: initPrimaryVariablesEvaluation() const { - for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) { - primary_variables_evaluation_[eqIdx] = - EvalWell::createVariable(numWellEq_ + numEq, primary_variables_[eqIdx], numEq + eqIdx); - } - } - - - - - - template - const typename StandardWell::EvalWell& - StandardWell:: - getBhp() const - { - return primary_variables_evaluation_[Bhp]; - } - - - - - - template - const typename StandardWell::EvalWell& - StandardWell:: - getWQTotal() const - { - return primary_variables_evaluation_[WQTotal]; - } - - - - - - template - typename StandardWell::EvalWell - StandardWell:: - getQs(const int comp_idx) const - { - // Note: currently, the WQTotal definition is still depends on Injector/Producer. - assert(comp_idx < num_components_); - - if (this->isInjector()) { // only single phase injection - double inj_frac = 0.0; - switch (this->wellEcl().injectorType()) { - case InjectorType::WATER: - if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) { - inj_frac = 1.0; - } - break; - case InjectorType::GAS: - if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent - inj_frac = wsolvent(); - } else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) { - inj_frac = has_solvent ? 1.0 - wsolvent() : 1.0; - } - break; - case InjectorType::OIL: - if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) { - inj_frac = 1.0; - } - break; - case InjectorType::MULTI: - // Not supported. - // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", - // "Multi phase injectors are not supported, requested for well " + name()); - break; - } - return inj_frac * primary_variables_evaluation_[WQTotal]; - } else { // producers - return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx); - } - } - - - - - - - template - typename StandardWell::EvalWell - StandardWell:: - wellVolumeFractionScaled(const int compIdx) const - { - - const int legacyCompIdx = ebosCompIdxToFlowCompIdx(compIdx); - const double scal = scalingFactor(legacyCompIdx); - if (scal > 0) - return wellVolumeFraction(compIdx) / scal; - - // the scaling factor may be zero for RESV controlled wells. - return wellVolumeFraction(compIdx); - } - - - - - - template - typename StandardWell::EvalWell - StandardWell:: - wellVolumeFraction(const unsigned compIdx) const - { - if (FluidSystem::numActivePhases() == 1) { - return EvalWell(numWellEq_ + numEq, 1.0); - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { - return primary_variables_evaluation_[WFrac]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { - return primary_variables_evaluation_[GFrac]; - } - - if (has_solvent && compIdx == (unsigned)contiSolventEqIdx) { - return primary_variables_evaluation_[SFrac]; - } - } - else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { - return primary_variables_evaluation_[GFrac]; - } - } - - // Oil or WATER fraction - EvalWell well_fraction(numWellEq_ + numEq, 1.0); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - well_fraction -= primary_variables_evaluation_[WFrac]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - well_fraction -= primary_variables_evaluation_[GFrac]; - } - - if (has_solvent) { - well_fraction -= primary_variables_evaluation_[SFrac]; - } - } - else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))) { - - well_fraction -= primary_variables_evaluation_[GFrac]; - } - - return well_fraction; - } - - - - - - template - typename StandardWell::EvalWell - StandardWell:: - wellSurfaceVolumeFraction(const int compIdx) const - { - - EvalWell sum_volume_fraction_scaled(numWellEq_ + numEq, 0.); - for (int idx = 0; idx < num_components_; ++idx) { - sum_volume_fraction_scaled += wellVolumeFractionScaled(idx); - } - - assert(sum_volume_fraction_scaled.value() != 0.); - - return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled; - } - - - - - - template - typename StandardWell::EvalWell - StandardWell:: - extendEval(const Eval& in) const - { - EvalWell out(numWellEq_ + numEq, in.value()); - for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { - out.setDerivative(eqIdx, in.derivative(eqIdx)); - } - return out; + this->StdWellEval::initPrimaryVariablesEvaluation(); } @@ -372,22 +114,20 @@ namespace Opm double& perf_vap_oil_rate, DeferredLogger& deferred_logger) const { - const auto& fs = intQuants.fluidState(); - const EvalWell pressure = extendEval(getPerfCellPressure(fs)); - const EvalWell rs = extendEval(fs.Rs()); - const EvalWell rv = extendEval(fs.Rv()); - std::vector b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0}); + const EvalWell pressure = this->extendEval(getPerfCellPressure(fs)); + const EvalWell rs = this->extendEval(fs.Rs()); + const EvalWell rv = this->extendEval(fs.Rv()); + std::vector b_perfcells_dense(num_components_, EvalWell{this->numWellEq_ + numEq, 0.0}); for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { continue; } - const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - b_perfcells_dense[compIdx] = extendEval(fs.invB(phaseIdx)); + b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx)); } if constexpr (has_solvent) { - b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); + b_perfcells_dense[contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor()); } if constexpr (has_zFraction) { @@ -397,142 +137,20 @@ namespace Opm b_perfcells_dense[gasCompIdx] += wsolvent()*intQuants.zPureInvFormationVolumeFactor().value(); } } - - // Pressure drawdown (also used to determine direction of flow) - const EvalWell well_pressure = bhp + this->perf_pressure_diffs_[perf]; - EvalWell drawdown = pressure - well_pressure; - - if constexpr (Base::has_polymermw) { - if (this->isInjector()) { - const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; - const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index]; - drawdown += skin_pressure; - } - } - - // producing perforations - if ( drawdown.value() > 0 ) { - //Do nothing if crossflow is not allowed - if (!allow_cf && this->isInjector()) { - return; - } - - // compute component volumetric rates at standard conditions - for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { - const EvalWell cq_p = - Tw * (mob[componentIdx] * drawdown); - cq_s[componentIdx] = b_perfcells_dense[componentIdx] * 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_sOil = cq_s[oilCompIdx]; - const EvalWell cq_sGas = cq_s[gasCompIdx]; - const EvalWell dis_gas = rs * cq_sOil; - const EvalWell vap_oil = rv * cq_sGas; - - cq_s[gasCompIdx] += dis_gas; - cq_s[oilCompIdx] += vap_oil; - - // recording the perforation solution gas rate and solution oil rates - if (this->isProducer()) { - perf_dis_gas_rate = dis_gas.value(); - perf_vap_oil_rate = vap_oil.value(); - } - } - - } else { - //Do nothing if crossflow is not allowed - if (!allow_cf && this->isProducer()) { - return; - } - - // Using total mobilities - EvalWell total_mob_dense = mob[0]; - for (int componentIdx = 1; componentIdx < num_components_; ++componentIdx) { - total_mob_dense += mob[componentIdx]; - } - - // injection perforations total volume rates - const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); - - // surface volume fraction of fluids within wellbore - std::vector cmix_s(num_components_, EvalWell{numWellEq_ + numEq}); - for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { - cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); - } - - // compute volume ratio between connection at standard conditions - EvalWell volumeRatio(numWellEq_ + numEq, 0.); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; - } - - if constexpr (has_solvent) { - volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx]; - } - - 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 - const EvalWell d = EvalWell(numWellEq_ + numEq, 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; - //std::cout << "tmp_oil " <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; - } - } - } + this->StdWellEval::computePerfRate(mob, + pressure, + bhp, + rs, + rv, + b_perfcells_dense, + Tw, + perf, + allow_cf, + has_polymermw, + cq_s, + perf_dis_gas_rate, + perf_vap_oil_rate, + deferred_logger); } @@ -587,9 +205,9 @@ namespace Opm auto& perf_rates = well_state.perfPhaseRates(this->index_of_well_); for (int perf = 0; perf < number_of_perforations_; ++perf) { // Calculate perforation quantities. - std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.0}); - EvalWell water_flux_s{numWellEq_ + numEq, 0.0}; - EvalWell cq_s_zfrac_effective{numWellEq_ + numEq, 0.0}; + std::vector cq_s(num_components_, {this->numWellEq_ + numEq, 0.0}); + EvalWell water_flux_s{this->numWellEq_ + numEq, 0.0}; + EvalWell cq_s_zfrac_effective{this->numWellEq_ + numEq, 0.0}; calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger); // Equation assembly for this perforation. @@ -609,7 +227,7 @@ namespace Opm this->resWell_[0][componentIdx] += cq_s_effective.value(); // assemble the jacobians - for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) { + for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+numEq); @@ -629,7 +247,7 @@ namespace Opm } if constexpr (has_zFraction) { - for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) { + for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { this->duneC_[0][cell_idx][pvIdx][contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+numEq); } } @@ -644,13 +262,13 @@ namespace Opm for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) { // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume // since all the rates are under surface condition - EvalWell resWell_loc(numWellEq_ + numEq, 0.0); + EvalWell resWell_loc(this->numWellEq_ + numEq, 0.0); if (FluidSystem::numActivePhases() > 1) { assert(dt > 0); - resWell_loc += (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; + resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt; } - resWell_loc -= getQs(componentIdx) * well_efficiency_factor_; - for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) { + resWell_loc -= this->getQs(componentIdx) * well_efficiency_factor_; + for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq); } this->resWell_[0][componentIdx] += resWell_loc.value(); @@ -658,7 +276,7 @@ namespace Opm const auto& summaryState = ebosSimulator.vanguard().summaryState(); const Schedule& schedule = ebosSimulator.vanguard().schedule(); - assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger); + this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger); // do the local inversion of D. @@ -685,10 +303,10 @@ namespace Opm DeferredLogger& deferred_logger) const { const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); - const EvalWell& bhp = getBhp(); + const EvalWell& bhp = this->getBhp(); const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); + std::vector mob(num_components_, {this->numWellEq_ + numEq, 0.}); getMobility(ebosSimulator, perf, mob, deferred_logger); double perf_dis_gas_rate = 0.; @@ -730,11 +348,11 @@ namespace Opm const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); // convert to reservoar conditions - EvalWell cq_r_thermal(numWellEq_ + numEq, 0.); + EvalWell cq_r_thermal(this->numWellEq_ + numEq, 0.); if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if(FluidSystem::waterPhaseIdx == phaseIdx) - cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); + cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx)); // remove dissolved gas and vapporized oil const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); @@ -742,16 +360,16 @@ namespace Opm // 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 - const EvalWell d = extendEval(1.0 - fs.Rv() * fs.Rs()); + const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs()); // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) if(FluidSystem::gasPhaseIdx == phaseIdx) - cq_r_thermal = (cq_s[gasCompIdx] - extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); + cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) ); // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) if(FluidSystem::oilPhaseIdx == phaseIdx) - cq_r_thermal = (cq_s[oilCompIdx] - extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); + cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) ); } else { - cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); + cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx)); } // change temperature for injecting fluids @@ -772,7 +390,7 @@ namespace Opm fs.setEnthalpy(phaseIdx, h); } // compute the thermal flux - cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx)); + cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx)); connectionRates[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal); } } @@ -784,7 +402,7 @@ namespace Opm if (this->isInjector()) { cq_s_poly *= wpolymer(); } else { - cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); + cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); } // Note. Efficiency factor is handled in the output layer auto& perf_rate_polymer = well_state.perfRatePolymer(this->index_of_well_); @@ -805,7 +423,7 @@ namespace Opm if (this->isInjector()) { cq_s_foam *= wfoam(); } else { - cq_s_foam *= extendEval(intQuants.foamConcentration()); + cq_s_foam *= this->extendEval(intQuants.foamConcentration()); } connectionRates[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam); } @@ -818,7 +436,7 @@ namespace Opm cq_s_zfrac_effective *= wsolvent(); } else if (cq_s_zfrac_effective.value() != 0.0) { const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value(); - cq_s_zfrac_effective *= extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume()); + cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume()); } auto& perf_rate_solvent = well_state.perfRateSolvent(this->index_of_well_); perf_rate_solvent[perf] = cq_s_zfrac_effective.value(); @@ -834,7 +452,7 @@ namespace Opm if (this->isInjector()) { cq_s_sm *= wsalt(); } else { - cq_s_sm *= extendEval(intQuants.fluidState().saltConcentration()); + cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration()); } // Note. Efficiency factor is handled in the output layer auto& perf_rate_brine = well_state.perfRateBrine(this->index_of_well_); @@ -852,68 +470,6 @@ namespace Opm - template - void - StandardWell::assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - DeferredLogger& deferred_logger) - { - EvalWell control_eq(numWellEq_ + numEq, 0.0); - - const auto& well = well_ecl_; - - auto getRates = [&]() { - std::vector rates(3, EvalWell(numWellEq_ + numEq, 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 injection rate. - const EvalWell injection_rate = getWQTotal(); - // 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. - const auto& inj_controls = well.injectionControls(summaryState); - 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. - const auto& prod_controls = well.productionControls(summaryState); - 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 - // TODO: we should use a different index system for the well equations - this->resWell_[0][Bhp] = control_eq.value(); - for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) { - this->invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + numEq); - } - } - - - - template void StandardWell:: @@ -939,10 +495,10 @@ 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()); + mob[contiSolventEqIdx] = this->extendEval(intQuants.solventMobility()); } } else { @@ -960,7 +516,7 @@ 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)); } // this may not work if viscosity and relperms has been modified? @@ -1013,54 +569,13 @@ namespace Opm const WellState& /* well_state */) const { const double dFLimit = param_.dwell_fraction_max_; - - const std::vector old_primary_variables = primary_variables_; - - // for injectors, very typical one of the fractions will be one, and it is easy to get zero value - // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here - const double relaxation_factor_fractions = (this->isProducer()) ? - relaxationFactorFractionsProducer(old_primary_variables, dwells) - : 1.0; - - // update the second and third well variable (The flux fractions) - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); - // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; - primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); - primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; - } - - if constexpr (has_solvent) { - const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; - const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); - primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; - } - - processFractions(); - - // updating the total rates Q_t - const double relaxation_factor_rate = this->relaxationFactorRate(old_primary_variables, dwells); - primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate; - - // updating the bottom hole pressure - { - const double dBHPLimit = param_.dbhp_max_rel_; - const int sign1 = dwells[0][Bhp] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit); - // 1e5 to make sure bhp will not be below 1bar - primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5); - } + const double dBHPLimit = param_.dbhp_max_rel_; + this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit); updateExtraPrimaryVariables(dwells); #ifndef NDEBUG - for (double v : primary_variables_) { + for (double v : this->primary_variables_) { assert(isfinite(v)); } #endif @@ -1078,118 +593,7 @@ namespace Opm { // for the water velocity and skin pressure if constexpr (Base::has_polymermw) { - if (this->isInjector()) { - for (int perf = 0; perf < number_of_perforations_; ++perf) { - const int wat_vel_index = Bhp + 1 + perf; - const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; - - const double relaxation_factor = 0.9; - const double dx_wat_vel = dwells[0][wat_vel_index]; - primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel; - - const double dx_pskin = dwells[0][pskin_index]; - primary_variables_[pskin_index] -= relaxation_factor * dx_pskin; - } - } - } - } - - - - - - template - void - StandardWell:: - processFractions() const - { - const auto pu = phaseUsage(); - std::vector F(number_of_phases_, 0.0); - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - F[pu.phase_pos[Oil]] = 1.0; - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - F[pu.phase_pos[Water]] = primary_variables_[WFrac]; - F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; - F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]]; - } - } - else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - F[pu.phase_pos[Water]] = 1.0; - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; - F[pu.phase_pos[Water]] -= F[pu.phase_pos[Gas]]; - } - } - else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - F[pu.phase_pos[Gas]] = 1.0; - } - - [[maybe_unused]] double F_solvent; - if constexpr (has_solvent) { - F_solvent = primary_variables_[SFrac]; - F[pu.phase_pos[Oil]] -= F_solvent; - } - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - if (F[Water] < 0.0) { - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]); - } - if constexpr (has_solvent) { - F_solvent /= (1.0 - F[pu.phase_pos[Water]]); - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]); - } - F[pu.phase_pos[Water]] = 0.0; - } - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - if (F[pu.phase_pos[Gas]] < 0.0) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); - } - if constexpr (has_solvent) { - F_solvent /= (1.0 - F[pu.phase_pos[Gas]]); - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]); - } - F[pu.phase_pos[Gas]] = 0.0; - } - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - if (F[pu.phase_pos[Oil]] < 0.0) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); - } - if constexpr (has_solvent) { - F_solvent /= (1.0 - F[pu.phase_pos[Oil]]); - } - F[pu.phase_pos[Oil]] = 0.0; - } - } - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - primary_variables_[WFrac] = F[pu.phase_pos[Water]]; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; - } - if constexpr (has_solvent) { - primary_variables_[SFrac] = F_solvent; + this->updatePrimaryVariablesPolyMW(dwells); } } @@ -1202,108 +606,11 @@ namespace Opm StandardWell:: updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const { - const PhaseUsage& pu = phaseUsage(); - std::vector F(number_of_phases_, 0.0); - [[maybe_unused]] double F_solvent = 0.0; - if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) { - const int oil_pos = pu.phase_pos[Oil]; - F[oil_pos] = 1.0; - - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; - F[water_pos] = primary_variables_[WFrac]; - F[oil_pos] -= F[water_pos]; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; - F[gas_pos] = primary_variables_[GFrac]; - F[oil_pos] -= F[gas_pos]; - } - - if constexpr (has_solvent) { - F_solvent = primary_variables_[SFrac]; - F[oil_pos] -= F_solvent; - } - } - else if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; - F[water_pos] = 1.0; - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; - F[gas_pos] = primary_variables_[GFrac]; - F[water_pos] -= F[gas_pos]; - } - } - else if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; - F[gas_pos] = 1.0; - } - - - // 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 scal = scalingFactor(p); - // for injection wells, there should only one non-zero scaling factor - if (scal > 0) { - F[p] /= scal ; - } else { - // this should only happens to injection wells - F[p] = 0.; - } - } - - // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. - // More testing is needed to make sure this is correct for well groups and THP. - if constexpr (has_solvent){ - F_solvent /= scalingFactor(contiSolventEqIdx); - F[pu.phase_pos[Gas]] += F_solvent; - } - - well_state.update_bhp(index_of_well_, primary_variables_[Bhp]); - - // calculate the phase rates based on the primary variables - // for producers, this is not a problem, while not sure for injectors here - if (this->isProducer()) { - const double g_total = primary_variables_[WQTotal]; - for (int p = 0; p < number_of_phases_; ++p) { - well_state.wellRates(index_of_well_)[p] = g_total * F[p]; - } - } else { // injectors - for (int p = 0; p < number_of_phases_; ++p) { - well_state.wellRates(index_of_well_)[p] = 0.0; - } - switch (this->wellEcl().injectorType()) { - case InjectorType::WATER: - well_state.wellRates(index_of_well_)[pu.phase_pos[Water]] = primary_variables_[WQTotal]; - break; - case InjectorType::GAS: - well_state.wellRates(index_of_well_)[pu.phase_pos[Gas]] = primary_variables_[WQTotal]; - break; - case InjectorType::OIL: - well_state.wellRates(index_of_well_)[pu.phase_pos[Oil]] = primary_variables_[WQTotal]; - break; - case InjectorType::MULTI: - // Not supported. - deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", - "Multi phase injectors are not supported, requested for well " + name()); - break; - } - } - - updateThp(well_state, deferred_logger); + this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger); // other primary variables related to polymer injectivity study if constexpr (Base::has_polymermw) { - if (this->isInjector()) { - auto& perf_water_velocity = well_state.perfWaterVelocity(this->index_of_well_); - auto& perf_skin_pressure = well_state.perfSkinPressure(this->index_of_well_); - for (int perf = 0; perf < number_of_perforations_; ++perf) { - perf_water_velocity[perf] = primary_variables_[Bhp + 1 + perf]; - perf_skin_pressure[perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf]; - } - } + this->updateWellStateFromPrimaryVariablesPolyMW(well_state); } } @@ -1311,55 +618,6 @@ namespace Opm - template - void - StandardWell:: - 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(well_state, rates, bhp, deferred_logger)); - - } - - - - - - template - void - StandardWell:: - updateWellStateWithTarget(const Simulator& ebos_simulator, - WellState& well_state, - DeferredLogger& deferred_logger) const - { - Base::updateWellStateWithTarget(ebos_simulator, well_state, deferred_logger); - } - - - - - template void StandardWell:: @@ -1378,7 +636,7 @@ namespace Opm std::fill(ipr_b_.begin(), ipr_b_.end(), 0.); for (int perf = 0; perf < number_of_perforations_; ++perf) { - std::vector mob(num_components_, {numWellEq_ + numEq, 0.0}); + std::vector mob(num_components_, {this->numWellEq_ + numEq, 0.0}); // TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration getMobility(ebos_simulator, perf, mob, deferred_logger); @@ -1561,7 +819,7 @@ namespace Opm const auto& fs = intQuants.fluidState(); const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value(); - const double bhp = getBhp().value(); + const double bhp = this->getBhp().value(); // Pressure drawdown (also used to determine direction of flow) const double well_pressure = bhp + this->perf_pressure_diffs_[perf]; @@ -1756,157 +1014,6 @@ namespace Opm - template - void - StandardWell:: - computeConnectionDensities(const std::vector& perfComponentRates, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf) - { - // Verify that we have consistent input. - const int nperf = number_of_perforations_; - const int num_comp = num_components_; - - // 1. Compute the flow (in surface volume units for each - // component) exiting up the wellbore from each perforation, - // taking into account flow from lower in the well, and - // in/out-flow at each perforation. - std::vector q_out_perf((nperf)*num_comp, 0.0); - - // Step 1 depends on the order of the perforations. Hence we need to - // do the modifications globally. - // Create and get the global perforation information and do this sequentially - // on each process - - const auto& factory = this->parallel_well_info_.getGlobalPerfContainerFactory(); - auto global_q_out_perf = factory.createGlobal(q_out_perf, num_comp); - auto global_perf_comp_rates = factory.createGlobal(perfComponentRates, num_comp); - - // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore - // Iterate over well perforations from bottom to top. - for (int perf = factory.numGlobalPerfs() - 1; perf >= 0; --perf) { - for (int component = 0; component < num_comp; ++component) { - auto index = perf * num_comp + component; - if (perf == factory.numGlobalPerfs() - 1) { - // This is the bottom perforation. No flow from below. - global_q_out_perf[index] = 0.0; - } else { - // Set equal to flow from below. - global_q_out_perf[index] = global_q_out_perf[index + num_comp]; - } - // Subtract outflow through perforation. - global_q_out_perf[index] -= global_perf_comp_rates[index]; - } - } - - // Copy the data back to local view - factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_comp); - - // 2. Compute the component mix at each perforation as the - // absolute values of the surface rates divided by their sum. - // Then compute volume ratios (formation factors) for each perforation. - // Finally compute densities for the segments associated with each perforation. - std::vector mix(num_comp,0.0); - std::vector x(num_comp); - std::vector surf_dens(num_comp); - - for (int perf = 0; perf < nperf; ++perf) { - // Find component mix. - const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf, - q_out_perf.begin() + num_comp*(perf+1), 0.0); - if (tot_surf_rate != 0.0) { - for (int component = 0; component < num_comp; ++component) { - mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); - } - } else if (num_comp == 1) { - mix[num_comp-1] = 1.0; - } else { - std::fill(mix.begin(), mix.end(), 0.0); - // No flow => use well specified fractions for mix. - if (this->isInjector()) { - switch (this->wellEcl().injectorType()) { - case InjectorType::WATER: - mix[FluidSystem::waterCompIdx] = 1.0; - break; - case InjectorType::GAS: - mix[FluidSystem::gasCompIdx] = 1.0; - break; - case InjectorType::OIL: - mix[FluidSystem::oilCompIdx] = 1.0; - break; - case InjectorType::MULTI: - // Not supported. - // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", - // "Multi phase injectors are not supported, requested for well " + name()); - break; - } - } else { - assert(this->isProducer()); - // For the frist perforation without flow we use the preferred phase to decide the mix initialization. - if (perf == 0) { // - switch (this->wellEcl().getPreferredPhase()) { - case Phase::OIL: - mix[FluidSystem::oilCompIdx] = 1.0; - break; - case Phase::GAS: - mix[FluidSystem::gasCompIdx] = 1.0; - break; - case Phase::WATER: - mix[FluidSystem::waterCompIdx] = 1.0; - break; - default: - // No others supported. - break; - } - // For the rest of the perforation without flow we use mix from the above perforation. - } else { - mix = x; - } - - } - } - // Compute volume ratio. - x = mix; - - // Subtract dissolved gas from oil phase and vapporized oil from gas phase - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - double rs = 0.0; - double rv = 0.0; - if (!rsmax_perf.empty() && mix[oilpos] > 1e-12) { - rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); - } - if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { - rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); - } - if (rs != 0.0) { - // Subtract gas in oil from gas mixture - x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv); - } - if (rv != 0.0) { - // Subtract oil in gas from oil mixture - x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);; - } - } - double volrat = 0.0; - for (int component = 0; component < num_comp; ++component) { - volrat += x[component] / b_perf[perf*num_comp+ component]; - } - for (int component = 0; component < num_comp; ++component) { - surf_dens[component] = surf_dens_perf[perf*num_comp+ component]; - } - - // Compute segment density. - this->perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; - } - } - - - - template ConvergenceReport StandardWell:: @@ -1922,43 +1029,13 @@ namespace Opm const double tol_wells = param_.tolerance_wells_; const double maxResidualAllowed = param_.max_residual_allowed_; - std::vector res(numWellEq_); - for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) { - // magnitude of the residual matters - res[eq_idx] = std::abs(this->resWell_[0][eq_idx]); - } - - std::vector well_flux_residual(num_components_); - - // Finish computation - for ( int compIdx = 0; compIdx < num_components_; ++compIdx ) - { - well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx]; - } - - ConvergenceReport report; - using CR = ConvergenceReport; - CR::WellFailure::Type type = CR::WellFailure::Type::MassBalance; - // checking if any NaN or too large residuals found - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - const int compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); - - if (std::isnan(well_flux_residual[compIdx])) { - report.setWellFailed({type, CR::Severity::NotANumber, compIdx, name()}); - } else if (well_flux_residual[compIdx] > maxResidualAllowed) { - report.setWellFailed({type, CR::Severity::TooLarge, compIdx, name()}); - } else if (well_flux_residual[compIdx] > tol_wells) { - report.setWellFailed({type, CR::Severity::Normal, compIdx, name()}); - } - } - - this->checkConvergenceControlEq(well_state, report, deferred_logger, maxResidualAllowed); - + std::vector res; + ConvergenceReport report = this->StdWellEval::getWellConvergence(well_state, + B_avg, + tol_wells, + maxResidualAllowed, + res, + deferred_logger); checkConvergenceExtraEqs(res, report); return report; @@ -2010,7 +1087,7 @@ namespace Opm return wellPICalc.connectionProdIndStandard(allPerfID, mobility); }; - std::vector mob(num_components_, {numWellEq_ + numEq, 0.0}); + std::vector mob(num_components_, {this->numWellEq_ + numEq, 0.0}); getMobility(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); const auto& fs = fluidState(subsetPerfID); @@ -2103,10 +1180,9 @@ namespace Opm } } - computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); this->computeConnectionPressureDelta(); - } @@ -2144,7 +1220,7 @@ namespace Opm // We assemble the well equations, then we check the convergence, // which is why we do not put the assembleWellEq here. BVectorWell dx_well(1); - dx_well[0].resize(numWellEq_); + dx_well[0].resize(this->numWellEq_); this->invDuneD_.mv(this->resWell_, dx_well); updateWellState(dx_well, well_state, deferred_logger); @@ -2164,25 +1240,11 @@ namespace Opm updatePrimaryVariables(well_state, deferred_logger); initPrimaryVariablesEvaluation(); computeWellConnectionPressures(ebosSimulator, well_state); - computeAccumWell(); + this->computeAccumWell(); } - template - void - StandardWell:: - computeAccumWell() - { - for (int eq_idx = 0; eq_idx < numWellConservationEq; ++eq_idx) { - F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value(); - } - } - - - - - template void StandardWell:: @@ -2229,58 +1291,6 @@ namespace Opm this->duneC_.mmtv(this->invDrw_, r); } -#if HAVE_CUDA || HAVE_OPENCL - template - void - StandardWell:: - addWellContribution(WellContributions& wellContribs) const - { - std::vector colIndices; - std::vector nnzValues; - colIndices.reserve(this->duneB_.nonzeroes()); - nnzValues.reserve(this->duneB_.nonzeroes()*numStaticWellEq * numEq); - - // duneC - for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC ) - { - colIndices.emplace_back(colC.index()); - for (int i = 0; i < numStaticWellEq; ++i) { - for (int j = 0; j < numEq; ++j) { - nnzValues.emplace_back((*colC)[i][j]); - } - } - } - wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->duneC_.nonzeroes()); - - // invDuneD - colIndices.clear(); - nnzValues.clear(); - colIndices.emplace_back(0); - for (int i = 0; i < numStaticWellEq; ++i) - { - for (int j = 0; j < numStaticWellEq; ++j) { - nnzValues.emplace_back(this->invDuneD_[0][0][i][j]); - } - } - wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1); - - // duneB - colIndices.clear(); - nnzValues.clear(); - for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB ) - { - colIndices.emplace_back(colB.index()); - for (int i = 0; i < numStaticWellEq; ++i) { - for (int j = 0; j < numEq; ++j) { - nnzValues.emplace_back((*colB)[i][j]); - } - } - } - wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->duneB_.nonzeroes()); - } -#endif - - template void StandardWell:: @@ -2309,7 +1319,7 @@ namespace Opm if (!this->isOperable() && !this->wellIsStopped()) return; BVectorWell xw(1); - xw[0].resize(numWellEq_); + xw[0].resize(this->numWellEq_); recoverSolutionWell(x, xw); updateWellState(xw, well_state, deferred_logger); @@ -2336,15 +1346,15 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); // flux for each perforation - std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); + std::vector mob(num_components_, {this->numWellEq_ + numEq, 0.}); getMobility(ebosSimulator, perf, mob, deferred_logger); double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); const double Tw = well_index_[perf] * trans_mult; - std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); + std::vector cq_s(num_components_, {this->numWellEq_ + numEq, 0.}); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRate(intQuants, mob, EvalWell(numWellEq_ + numEq, bhp), Tw, perf, allow_cf, + computePerfRate(intQuants, mob, EvalWell(this->numWellEq_ + numEq, bhp), Tw, perf, allow_cf, cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); for(int p = 0; p < np; ++p) { @@ -2469,10 +1479,12 @@ namespace Opm std::vector &potentials, double alq) const { - /*double bhp =*/ computeWellRatesAndBhpWithThpAlqProd( - ebos_simulator, summary_state, - deferred_logger, potentials, alq - ); + /*double bhp =*/ + computeWellRatesAndBhpWithThpAlqProd(ebos_simulator, + summary_state, + deferred_logger, + potentials, + alq); } template @@ -2578,109 +1590,22 @@ namespace Opm StandardWell:: updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const { + this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger); if (!this->isOperable() && !this->wellIsStopped()) return; - const int well_index = index_of_well_; - const int np = number_of_phases_; - const auto& pu = phaseUsage(); - - // the weighted total well rate - double total_well_rate = 0.0; - for (int p = 0; p < np; ++p) { - total_well_rate += scalingFactor(p) * well_state.wellRates(well_index)[p]; - } - - // Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate - // under surface condition is used here - if (this->isInjector()) { - switch (this->wellEcl().injectorType()) { - case InjectorType::WATER: - primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Water]]; - break; - case InjectorType::GAS: - primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Gas]]; - break; - case InjectorType::OIL: - primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Oil]]; - break; - case InjectorType::MULTI: - // Not supported. - deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", - "Multi phase injectors are not supported, requested for well " + name()); - break; - } - } else { - primary_variables_[WQTotal] = total_well_rate; - } - - if (std::abs(total_well_rate) > 0.) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates(well_index)[pu.phase_pos[Water]] / total_well_rate; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates(well_index)[pu.phase_pos[Gas]] - - (has_solvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ; - } - if constexpr (has_solvent) { - primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; - } - } else { // total_well_rate == 0 - if (this->isInjector()) { - auto phase = well_ecl_.getInjectionProperties().injectorType; - // only single phase injection handled - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - if (phase == InjectorType::WATER) { - primary_variables_[WFrac] = 1.0; - } else { - primary_variables_[WFrac] = 0.0; - } - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - if (phase == InjectorType::GAS) { - primary_variables_[GFrac] = 1.0; - if constexpr (has_solvent) { - primary_variables_[GFrac] = 1.0 - wsolvent(); - primary_variables_[SFrac] = wsolvent(); - } - } else { - primary_variables_[GFrac] = 0.0; - } - } - - // TODO: it is possible to leave injector as a oil well, - // when F_w and F_g both equals to zero, not sure under what kind of circumstance - // this will happen. - } else if (this->isProducer()) { // producers - // TODO: the following are not addressed for the solvent case yet - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - primary_variables_[WFrac] = 1.0 / np; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - primary_variables_[GFrac] = 1.0 / np; - } - } else { - OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger); - } - } - - - // BHP - primary_variables_[Bhp] = well_state.bhp(index_of_well_); - // other primary variables related to polymer injection if constexpr (Base::has_polymermw) { if (this->isInjector()) { const auto& water_velocity = well_state.perfWaterVelocity(this->index_of_well_); const auto& skin_pressure = well_state.perfSkinPressure(this->index_of_well_); for (int perf = 0; perf < number_of_perforations_; ++perf) { - primary_variables_[Bhp + 1 + perf] = water_velocity[perf]; - primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = skin_pressure[perf]; + this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf]; + this->primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = skin_pressure[perf]; } } } #ifndef NDEBUG - for (double v : primary_variables_) { + for (double v : this->primary_variables_) { assert(isfinite(v)); } #endif @@ -2710,7 +1635,7 @@ namespace Opm { const int cell_idx = well_cells_[perf]; const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - const EvalWell polymer_concentration = extendEval(int_quant.polymerConcentration()); + const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration()); // TODO: not sure should based on the well type or injecting/producing peforations // it can be different for crossflow @@ -2718,7 +1643,7 @@ namespace Opm // assume fully mixing within injecting wellbore const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex()); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - mob[waterCompIdx] /= (extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) ); + mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) ); } if (PolymerModule::hasPlyshlog()) { @@ -2730,9 +1655,9 @@ namespace Opm // compute the well water velocity with out shear effects. // TODO: do we need to turn on crossflow here? const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator); - const EvalWell& bhp = getBhp(); + const EvalWell& bhp = this->getBhp(); - std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); + std::vector cq_s(num_components_, {this->numWellEq_ + numEq, 0.}); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier(int_quant, cell_idx); @@ -2745,12 +1670,12 @@ namespace Opm const auto& scaled_drainage_info = material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx); const double swcr = scaled_drainage_info.Swcr; - const EvalWell poro = extendEval(int_quant.porosity()); - const EvalWell sw = extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx)); + const EvalWell poro = this->extendEval(int_quant.porosity()); + const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx)); // guard against zero porosity and no water const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - EvalWell water_velocity = cq_s[waterCompIdx] / denom * extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx)); + EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx)); if (PolymerModule::hasShrate()) { // the equation for the water velocity conversion for the wells and reservoir are from different version @@ -2793,53 +1718,6 @@ namespace Opm - template - double - StandardWell:: - relaxationFactorFractionsProducer(const std::vector& primary_variables, - const BVectorWell& dwells) - { - // TODO: not considering solvent yet - // 0.95 is a experimental value, which remains to be optimized - double relaxation_factor = 1.0; - - if (FluidSystem::numActivePhases() > 1) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const double relaxation_factor_w = StandardWellGeneric:: - relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const double relaxation_factor_g = StandardWellGeneric:: - relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); - } - - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will - // not be negative oil fraction later - const double original_sum = primary_variables[WFrac] + primary_variables[GFrac]; - const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor; - const double possible_updated_sum = original_sum - relaxed_update; - - if (possible_updated_sum > 1.0) { - assert(relaxed_update != 0.); - - const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95; - relaxation_factor *= further_relaxation_factor; - } - } - - assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0); - } - return relaxation_factor; - } - - - - - template typename StandardWell::EvalWell StandardWell:: @@ -2853,9 +1731,9 @@ namespace Opm OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger); } const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id); - const EvalWell throughput_eval(numWellEq_ + numEq, throughput); + const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); // the skin pressure when injecting water, which also means the polymer concentration is zero - EvalWell pskin_water(numWellEq_ + numEq, 0.0); + EvalWell pskin_water(this->numWellEq_ + numEq, 0.0); pskin_water = water_table_func.eval(throughput_eval, water_velocity); return pskin_water; } else { @@ -2888,9 +1766,9 @@ namespace Opm } const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id); const double reference_concentration = skprpolytable.refConcentration; - const EvalWell throughput_eval(numWellEq_ + numEq, throughput); + const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); // the skin pressure when injecting water, which also means the polymer concentration is zero - EvalWell pskin_poly(numWellEq_ + numEq, 0.0); + EvalWell pskin_poly(this->numWellEq_ + numEq, 0.0); pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs); if (poly_inj_conc == reference_concentration) { return sign * pskin_poly; @@ -2919,8 +1797,8 @@ namespace Opm if constexpr (Base::has_polymermw) { const int table_id = well_ecl_.getPolymerProperties().m_plymwinjtable; const auto& table_func = PolymerModule::getPlymwinjTable(table_id); - const EvalWell throughput_eval(numWellEq_ + numEq, throughput); - EvalWell molecular_weight(numWellEq_ + numEq, 0.); + const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); + EvalWell molecular_weight(this->numWellEq_ + numEq, 0.); if (wpolymer() == 0.) { // not injecting polymer return molecular_weight; } @@ -2945,7 +1823,7 @@ namespace Opm if (this->isInjector()) { auto& perf_water_throughput = well_state.perfThroughput(this->index_of_well_); for (int perf = 0; perf < number_of_perforations_; ++perf) { - const double perf_water_vel = primary_variables_[Bhp + 1 + perf]; + const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf]; // we do not consider the formation damage due to water flowing from reservoir into wellbore if (perf_water_vel > 0.) { perf_water_throughput[perf] += perf_water_vel * dt; @@ -2969,14 +1847,14 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quants.fluidState(); - const EvalWell b_w = extendEval(fs.invB(FluidSystem::waterPhaseIdx)); + const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx)); const double area = M_PI * bore_diameters_[perf] * perf_length_[perf]; const int wat_vel_index = Bhp + 1 + perf; const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); // water rate is update to use the form from water velocity, since water velocity is // a primary variable now - cq_s[water_comp_idx] = area * primary_variables_evaluation_[wat_vel_index] * b_w; + cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w; } @@ -2994,29 +1872,29 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quants.fluidState(); - const EvalWell b_w = extendEval(fs.invB(FluidSystem::waterPhaseIdx)); + const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx)); const EvalWell water_flux_r = water_flux_s / b_w; const double area = M_PI * bore_diameters_[perf] * perf_length_[perf]; const EvalWell water_velocity = water_flux_r / area; const int wat_vel_index = Bhp + 1 + perf; // equation for the water velocity - const EvalWell eq_wat_vel = primary_variables_evaluation_[wat_vel_index] - water_velocity; + const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity; this->resWell_[0][wat_vel_index] = eq_wat_vel.value(); const auto& perf_water_throughput = well_state.perfThroughput(this->index_of_well_); const double throughput = perf_water_throughput[perf]; const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; - EvalWell poly_conc(numWellEq_ + numEq, 0.0); + EvalWell poly_conc(this->numWellEq_ + numEq, 0.0); poly_conc.setValue(wpolymer()); // equation for the skin pressure - const EvalWell eq_pskin = primary_variables_evaluation_[pskin_index] - - pskin(throughput, primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger); + const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index] + - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger); this->resWell_[0][pskin_index] = eq_pskin.value(); - for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) { + for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+numEq); this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+numEq); } @@ -3063,7 +1941,7 @@ namespace Opm EvalWell cq_s_polymw = cq_s_poly; if (this->isInjector()) { const int wat_vel_index = Bhp + 1 + perf; - const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index]; + const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index]; if (water_velocity > 0.) { // injecting const auto& perf_water_throughput = well_state.perfThroughput(this->index_of_well_); const double throughput = perf_water_throughput[perf]; @@ -3076,7 +1954,7 @@ namespace Opm } } else if (this->isProducer()) { if (cq_s_polymw < 0.) { - cq_s_polymw *= extendEval(int_quants.polymerMoleWeight() ); + cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() ); } else { // we do not consider the molecular weight from the polymer // re-injecting back through producer @@ -3099,8 +1977,10 @@ namespace Opm const SummaryState& summary_state, DeferredLogger& deferred_logger) const { - return computeBhpAtThpLimitProdWithAlq( - ebos_simulator, summary_state, deferred_logger, getALQ(well_state)); + return computeBhpAtThpLimitProdWithAlq(ebos_simulator, + summary_state, + deferred_logger, + getALQ(well_state)); } template @@ -3205,15 +2085,15 @@ namespace Opm DeferredLogger& deferred_logger) const { // Calculate the rates that follow from the current primary variables. - std::vector well_q_s(num_components_, {numWellEq_ + numEq, 0.}); - const EvalWell& bhp = getBhp(); + std::vector well_q_s(num_components_, {this->numWellEq_ + numEq, 0.}); + const EvalWell& bhp = this->getBhp(); const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); for (int perf = 0; perf < number_of_perforations_; ++perf) { const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); + std::vector mob(num_components_, {this->numWellEq_ + numEq, 0.}); getMobility(ebosSimulator, perf, mob, deferred_logger); - std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); + std::vector cq_s(num_components_, {this->numWellEq_ + numEq, 0.}); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); diff --git a/opm/simulators/wells/WellInterfaceEval.cpp b/opm/simulators/wells/WellInterfaceEval.cpp index c762329a0..b762e6579 100644 --- a/opm/simulators/wells/WellInterfaceEval.cpp +++ b/opm/simulators/wells/WellInterfaceEval.cpp @@ -66,7 +66,7 @@ getGroupInjectionControl(const Group& group, const EvalWell& injection_rate, EvalWell& control_eq, double efficiencyFactor, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger) const { // Setting some defaults to silence warnings below. // Will be overwritten in the switch statement. @@ -183,7 +183,7 @@ getGroupProductionControl(const Group& group, const EvalWell& bhp, const std::vector& rates, EvalWell& control_eq, - double efficiencyFactor) + double efficiencyFactor) const { const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); if (currentGroupControl == Group::ProductionCMode::FLD || @@ -278,7 +278,7 @@ assembleControlEqProd_(const WellState& well_state, const std::vector& rates, // Always 3 canonical rates. const std::function& bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger) const { auto current = well_state.currentProductionControl(baseif_.indexOfWell()); const auto& pu = baseif_.phaseUsage(); @@ -392,7 +392,7 @@ assembleControlEqInj_(const WellState& well_state, const EvalWell& injection_rate, const std::function& bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger) const { auto current = well_state.currentInjectionControl(baseif_.indexOfWell()); const InjectorType injectorType = controls.injector_type; @@ -514,7 +514,7 @@ assembleControlEqProd_<__VA_ARGS__>(const WellState&, \ const std::vector<__VA_ARGS__>&, \ const std::function<__VA_ARGS__()>&, \ __VA_ARGS__&, \ - DeferredLogger&); \ + DeferredLogger&) const; \ template void WellInterfaceEval:: \ assembleControlEqInj_<__VA_ARGS__>(const WellState&, \ const GroupState&, \ @@ -525,7 +525,7 @@ assembleControlEqInj_<__VA_ARGS__>(const WellState&, \ const __VA_ARGS__&, \ const std::function<__VA_ARGS__()>&, \ __VA_ARGS__&, \ - DeferredLogger&); \ + DeferredLogger&) const; \ template __VA_ARGS__ WellInterfaceEval:: \ calculateBhpFromThp<__VA_ARGS__>(const WellState&, \ const std::vector<__VA_ARGS__>&, \ diff --git a/opm/simulators/wells/WellInterfaceEval.hpp b/opm/simulators/wells/WellInterfaceEval.hpp index d36f04a64..1a558a1d2 100644 --- a/opm/simulators/wells/WellInterfaceEval.hpp +++ b/opm/simulators/wells/WellInterfaceEval.hpp @@ -48,33 +48,14 @@ class WellInterfaceEval { static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Gas = BlackoilPhases::Vapour; -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); - - - 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); +public: + template + EvalWell calculateBhpFromThp(const WellState& well_state, + const std::vector& rates, + const Well& well, + const SummaryState& summaryState, + const double rho, + DeferredLogger& deferred_logger) const; template void assembleControlEqProd(const WellState& well_state, @@ -86,7 +67,7 @@ protected: const std::vector& rates, // Always 3 canonical rates. BhpFromThpFunc bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger) const { std::function eval = [&bhp_from_thp]() { return bhp_from_thp(); }; assembleControlEqProd_(well_state, @@ -111,7 +92,7 @@ protected: const std::vector& rates, // Always 3 canonical rates. const std::function& bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger); + DeferredLogger& deferred_logger) const; template void assembleControlEqInj(const WellState& well_state, @@ -123,7 +104,7 @@ protected: const EvalWell& injection_rate, BhpFromThpFunc bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger) const { std::function eval = [&bhp_from_thp]() { return bhp_from_thp(); }; assembleControlEqInj_(well_state, @@ -148,15 +129,59 @@ protected: const EvalWell& injection_rate, const std::function& bhp_from_thp, EvalWell& control_eq, - DeferredLogger& deferred_logger); + DeferredLogger& deferred_logger) const; - template - EvalWell calculateBhpFromThp(const WellState& well_state, - const std::vector& rates, - const Well& well, - const SummaryState& summaryState, - const double rho, - DeferredLogger& deferred_logger) const; +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 2e55cc2d5..b7aae6123 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -118,6 +118,10 @@ public: return guide_rate_; } + int numComponents() const{ + return num_components_; + } + int numPhases() const { return number_of_phases_; } @@ -156,13 +160,13 @@ public: double getTHPConstraint(const SummaryState& summaryState) const; double getALQ(const WellState& well_state) const; + double wsolvent() const; -protected: // whether a well is specified with a non-zero and valid VFP table number bool isVFPActive(DeferredLogger& deferred_logger) const; +protected: bool getAllowCrossFlow() const; - double wsolvent() const; double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; void updateWellTestStatePhysical(const WellState& well_state, const double simulation_time, diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index bf44c47a1..284321878 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -129,6 +129,7 @@ INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,fa INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>) INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u>) // Blackoil INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) @@ -138,6 +139,7 @@ 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>) // Alternative indices INSTANCE(EclAlternativeBlackOilIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) diff --git a/opm/simulators/wells/WellInterfaceIndices.hpp b/opm/simulators/wells/WellInterfaceIndices.hpp index b024a7c7a..98e85405a 100644 --- a/opm/simulators/wells/WellInterfaceIndices.hpp +++ b/opm/simulators/wells/WellInterfaceIndices.hpp @@ -38,6 +38,10 @@ public: using WellInterfaceFluidSystem::Oil; using WellInterfaceFluidSystem::Water; + int flowPhaseToEbosCompIdx(const int phaseIdx) const; + int ebosCompIdxToFlowCompIdx(const unsigned compIdx) const; + double scalingFactor(const int phaseIdx) const; + protected: WellInterfaceIndices(const Well& well, const ParallelWellInfo& parallel_well_info, @@ -48,11 +52,6 @@ protected: const int num_phases, const int index_of_well, const std::vector& perf_data); - - int flowPhaseToEbosCompIdx( const int phaseIdx ) const; - int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const; - - double scalingFactor(const int phaseIdx) const; }; }