diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0ccd767c1..af9783ff1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/GlobalWellInfo.cpp opm/simulators/wells/GroupState.cpp opm/simulators/wells/MSWellHelpers.cpp + opm/simulators/wells/MultisegmentWellAssemble.cpp opm/simulators/wells/MultisegmentWellEquations.cpp opm/simulators/wells/MultisegmentWellEval.cpp opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -372,6 +373,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/MSWellHelpers.hpp opm/simulators/wells/MultisegmentWell.hpp opm/simulators/wells/MultisegmentWell_impl.hpp + opm/simulators/wells/MultisegmentWellAssemble.hpp opm/simulators/wells/MultisegmentWellEquations.hpp opm/simulators/wells/MultisegmentWellEval.hpp opm/simulators/wells/MultisegmentWellGeneric.hpp diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp new file mode 100644 index 000000000..5a8382bfd --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -0,0 +1,370 @@ +/* + 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 + +namespace Opm { + +//! \brief Class administering assembler access to equation system. +template +class MultisegmentWellEquationAccess { +public: + //! \brief Constructor initializes reference to the equation system. + MultisegmentWellEquationAccess(MultisegmentWellEquations& eqns) + : eqns_(eqns) + {} + + using BVectorWell = typename MultisegmentWellEquations::BVectorWell; + using DiagMatWell = typename MultisegmentWellEquations::DiagMatWell; + using OffDiatMatWell = typename MultisegmentWellEquations::OffDiagMatWell; + + //! \brief Returns a reference to residual vector. + BVectorWell& residual() + { + return eqns_.resWell_; + } + + //! \brief Returns a reference to B matrix. + OffDiatMatWell& B() + { + return eqns_.duneB_; + } + + //! \brief Returns a reference to C matrix. + OffDiatMatWell& C() + { + return eqns_.duneC_; + } + + //! \brief Returns a reference to D matrix. + DiagMatWell& D() + { + return eqns_.duneD_; + } + +private: + MultisegmentWellEquations& eqns_; //!< Reference to equation system +}; + +template +void MultisegmentWellAssemble:: +assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + const double rho, + const EvalWell& wqTotal, + const EvalWell& bhp, + const std::function& getQs, + Equations& eqns1, + DeferredLogger& deferred_logger) const +{ + static constexpr int Gas = BlackoilPhases::Vapour; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Water = BlackoilPhases::Aqua; + + EvalWell control_eq(0.0); + + const auto& well = well_.wellEcl(); + + auto getRates = [&]() { + std::vector rates(3, 0.0); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)); + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); + } + return rates; + }; + + if (well_.wellIsStopped()) { + control_eq = wqTotal; + } else if (well_.isInjector() ) { + // Find scaling factor to get injection rate, + const InjectorType injectorType = inj_controls.injector_type; + double scaling = 1.0; + const auto& pu = well_.phaseUsage(); + switch (injectorType) { + case InjectorType::WATER: + { + scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]); + break; + } + case InjectorType::OIL: + { + scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]); + break; + } + case InjectorType::GAS: + { + scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]); + break; + } + default: + throw("Expected WATER, OIL or GAS as type for injectors " + well.name()); + } + const EvalWell injection_rate = wqTotal / scaling; + // Setup function for evaluation of BHP from THP (used only if needed). + std::function bhp_from_thp = [&]() { + const auto rates = getRates(); + return WellBhpThpCalculator(well_).calculateBhpFromThp(well_state, + rates, + well, + summaryState, + rho, + deferred_logger); + }; + // Call generic implementation. + WellAssemble(well_).assembleControlEqInj(well_state, + group_state, + schedule, + summaryState, + inj_controls, + bhp, + 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). + std::function bhp_from_thp = [&]() { + return WellBhpThpCalculator(well_).calculateBhpFromThp(well_state, + rates, + well, + summaryState, + rho, + deferred_logger); + }; + // Call generic implementation. + WellAssemble(well_).assembleControlEqProd(well_state, + group_state, + schedule, + summaryState, + prod_controls, + bhp, + rates, + bhp_from_thp, + control_eq, + deferred_logger); + } + + MultisegmentWellEquationAccess eqns(eqns1); + // using control_eq to update the matrix and residuals + eqns.residual()[0][SPres] = control_eq.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.D()[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); + } +} + +template +void MultisegmentWellAssemble:: +assemblePressureLoss(const int seg, + const int seg_upwind, + const EvalWell& accelerationPressureLoss, + Equations& eqns1) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][SPres] -= accelerationPressureLoss.value(); + eqns.D()[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); + eqns.D()[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq); + if constexpr (has_wfrac_variable) { + eqns.D()[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); + } + if constexpr (has_gfrac_variable) { + eqns.D()[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); + } +} + +template +void MultisegmentWellAssemble:: +assemblePressureEq(const int seg, + const int seg_upwind, + const int outlet_segment_index, + const EvalWell& pressure_equation, + const EvalWell& outlet_pressure, + Equations& eqns1, + bool wfrac, + bool gfrac) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][SPres] = pressure_equation.value(); + eqns.D()[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); + eqns.D()[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); + if (wfrac) { + eqns.D()[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); + } + if (gfrac) { + eqns.D()[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); + } + + // contribution from the outlet segment + eqns.residual()[seg][SPres] -= outlet_pressure.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.D()[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); + } +} + +template +void MultisegmentWellAssemble:: +assembleTrivialEq(const int seg, + const Scalar value, + Equations& eqns1) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][SPres] = value; + eqns.D()[seg][seg][SPres][WQTotal] = 1.; +} + +template +void MultisegmentWellAssemble:: +assembleAccumulationTerm(const int seg, + const int comp_idx, + const EvalWell& accumulation_term, + Equations& eqns1) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][comp_idx] += accumulation_term.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.D()[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq); + } +} + +template +void MultisegmentWellAssemble:: +assembleOutflowTerm(const int seg, + const int seg_upwind, + const int comp_idx, + const EvalWell& segment_rate, + Equations& eqns1) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][comp_idx] -= segment_rate.value(); + eqns.D()[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + eqns.D()[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + eqns.D()[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); + } + // pressure derivative should be zero +} + +template +void MultisegmentWellAssemble:: +assembleInflowTerm(const int seg, + const int inlet, + const int inlet_upwind, + const int comp_idx, + const EvalWell& inlet_rate, + Equations& eqns1) const + { + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][comp_idx] += inlet_rate.value(); + eqns.D()[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + eqns.D()[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + eqns.D()[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); + } + // pressure derivative should be zero +} + +template +void MultisegmentWellAssemble:: +assemblePerforationEq(const int seg, + const int cell_idx, + const int comp_idx, + const EvalWell& cq_s_effective, + Equations& eqns1) const +{ + MultisegmentWellEquationAccess eqns(eqns1); + // subtract sum of phase fluxes in the well equations. + eqns.residual()[seg][comp_idx] += cq_s_effective.value(); + + // assemble the jacobians + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + // also need to consider the efficiency factor when manipulating the jacobians. + eqns.C()[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // input in transformed matrix + + // the index name for the D should be eq_idx / pv_idx + eqns.D()[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq); + } + + for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) { + // also need to consider the efficiency factor when manipulating the jacobians. + eqns.B()[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); + } +} + +#define INSTANCE(...) \ +template class MultisegmentWellAssemble,__VA_ARGS__,double>; + +// One phase +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) + +// Two phase +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) + +// Blackoil +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,1u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>) +INSTANCE(BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>) + +} diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp new file mode 100644 index 000000000..693aa0c0e --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -0,0 +1,140 @@ +/* + 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_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED + +#include +#include + +#include + +namespace Opm +{ + +class DeferredLogger; +class GroupState; +template class MultisegmentWellEquations; +class Schedule; +class SummaryState; +template class WellInterfaceIndices; +class WellState; + +//! \brief Class handling assemble of the equation system for MultisegmentWell. +template +class MultisegmentWellAssemble +{ + static constexpr bool has_water = (Indices::waterSwitchIdx >= 0); + static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0); + static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0; + + // In the implementation, one should use has_wfrac_variable + // rather than has_water to check if you should do something + // with the variable at the WFrac location, similar for GFrac. + static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1; + static constexpr bool has_gfrac_variable = has_gas && has_oil; + + static constexpr int WQTotal = 0; + static constexpr int WFrac = has_wfrac_variable ? 1 : -1000; + static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000; + static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1; + +public: + static constexpr int numWellEq = Indices::numPhases+1; + using Equations = MultisegmentWellEquations; + using EvalWell = DenseAd::Evaluation; + //! \brief Constructor initializes reference to well. + MultisegmentWellAssemble(const WellInterfaceIndices& well) + : well_(well) + {} + + //! \brief Assemble control equation. + void assembleControlEq(const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + const double rho, + const EvalWell& wqTotal, + const EvalWell& bhp, + const std::function& getQs, + Equations& eqns, + DeferredLogger& deferred_logger) const; + + + //! \brief Assemble pressure loss term. + void assemblePressureLoss(const int seg, + const int seg_upwind, + const EvalWell& accelerationPressureLoss, + Equations& eqns) const; + + //! \brief Assemble pressure terms. + void assemblePressureEq(const int seg, + const int seg_upwind, + const int outlet_segment_index, + const EvalWell& pressure_equation, + const EvalWell& outlet_pressure, + Equations& eqns, + bool wfrac = has_wfrac_variable, + bool gfrac = has_gfrac_variable) const; + + //! \brief Assembles a trivial equation. + void assembleTrivialEq(const int seg, + const Scalar value, + Equations& eqns) const; + + //! \brief Assemble accumulation term. + void assembleAccumulationTerm(const int seg, + const int comp_idx, + const EvalWell& accumulation_term, + Equations& eqns1) const; + + //! \brief Assemble outflow term. + void assembleOutflowTerm(const int seg, + const int seg_upwind, + const int comp_idx, + const EvalWell& segment_rate, + Equations& eqns1) const; + + //! \brief Assemble inflow term. + void assembleInflowTerm(const int seg, + const int inlet, + const int inlet_upwind, + const int comp_idx, + const EvalWell& inlet_rate, + Equations& eqns) const; + + //! \brief Assemble equation for a perforation. + void assemblePerforationEq(const int seg, + const int cell_idx, + const int comp_idx, + const EvalWell& cq_s_effective, + Equations& eqns) const; + +private: + const WellInterfaceIndices& well_; //!< Reference to well +}; + +} + +#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index 9bdb79d0a..671d61459 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -36,6 +36,7 @@ template class UMFPack; namespace Opm { +template class MultisegmentWellEquationAccess; template class MultisegmentWellGeneric; class WellContributions; class WellInterfaceGeneric; @@ -110,6 +111,14 @@ public: const int seg_pressure_var_ind, const WellState& well_state) const; + //! \brief Returns a const reference to the residual. + const BVectorWell& residual() const + { + return resWell_; + } + + private: + friend class MultisegmentWellEquationAccess; // two off-diagonal matrices OffDiagMatWell duneB_; OffDiagMatWell duneC_; @@ -119,12 +128,11 @@ public: /// \brief solver for diagonal matrix /// /// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials... - mutable std::shared_ptr > duneDSolver_; + mutable std::shared_ptr> duneDSolver_; // residuals of the well equations BVectorWell resWell_; -private: const MultisegmentWellGeneric& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 2514ef879..3db453e69 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -112,7 +113,7 @@ getWellConvergence(const WellState& well_state, std::vector> abs_residual(this->numberOfSegments(), std::vector(numWellEq, 0.0)); for (int seg = 0; seg < this->numberOfSegments(); ++seg) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - abs_residual[seg][eq_idx] = std::abs(linSys_.resWell_[seg][eq_idx]); + abs_residual[seg][eq_idx] = std::abs(linSys_.residual()[seg][eq_idx]); } } @@ -177,7 +178,7 @@ getWellConvergence(const WellState& well_state, tolerance_wells, tolerance_wells, max_residual_allowed}, - std::abs(linSys_.resWell_[0][SPres]), + std::abs(linSys_.residual()[0][SPres]), report, deferred_logger); @@ -1107,120 +1108,6 @@ getSegmentSurfaceVolume(const EvalWell& temperature, return volume / vol_ratio; } -template -void -MultisegmentWellEval:: -assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::InjectionControls& inj_controls, - const Well::ProductionControls& prod_controls, - const double rho, - DeferredLogger& deferred_logger) -{ - static constexpr int Gas = BlackoilPhases::Vapour; - static constexpr int Oil = BlackoilPhases::Liquid; - static constexpr int Water = BlackoilPhases::Aqua; - - EvalWell control_eq(0.0); - - const auto& well = baseif_.wellEcl(); - - auto getRates = [&]() { - std::vector rates(3, 0.0); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)); - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx)); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); - } - return rates; - }; - - if (baseif_.wellIsStopped()) { - control_eq = getWQTotal(); - } else if (baseif_.isInjector() ) { - // Find scaling factor to get injection rate, - const InjectorType injectorType = inj_controls.injector_type; - double scaling = 1.0; - const auto& pu = baseif_.phaseUsage(); - switch (injectorType) { - case InjectorType::WATER: - { - scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]); - break; - } - case InjectorType::OIL: - { - scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]); - break; - } - case InjectorType::GAS: - { - scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]); - break; - } - default: - throw("Expected WATER, OIL or GAS as type for injectors " + well.name()); - } - const EvalWell injection_rate = getWQTotal() / scaling; - // Setup function for evaluation of BHP from THP (used only if needed). - std::function bhp_from_thp = [&]() { - const auto rates = getRates(); - return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state, - rates, - well, - summaryState, - rho, - deferred_logger); - }; - // Call generic implementation. - WellAssemble(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). - std::function bhp_from_thp = [&]() { - return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state, - rates, - well, - summaryState, - rho, - deferred_logger); - }; - // Call generic implementation. - WellAssemble(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 - linSys_.resWell_[0][SPres] = control_eq.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - linSys_.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); - } -} - template void MultisegmentWellEval:: @@ -1262,15 +1149,8 @@ handleAccelerationPressureLoss(const int seg, auto& segments = well_state.well(baseif_.indexOfWell()).segments; segments.pressure_drop_accel[seg] = accelerationPressureLoss.value(); - linSys_.resWell_[seg][SPres] -= accelerationPressureLoss.value(); - linSys_.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); - linSys_.duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq); - if (has_wfrac_variable) { - linSys_.duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); - } - if (has_gfrac_variable) { - linSys_.duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); - } + MultisegmentWellAssemble(baseif_). + assemblePressureLoss(seg, seg_upwind, accelerationPressureLoss, linSys_); } template @@ -1300,25 +1180,14 @@ assembleDefaultPressureEq(const int seg, segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); } - linSys_.resWell_[seg][SPres] = pressure_equation.value(); - const int seg_upwind = upwinding_segments_[seg]; - linSys_.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); - linSys_.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); - if (has_wfrac_variable) { - linSys_.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); - } - if (has_gfrac_variable) { - linSys_.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); - } - // contribution from the outlet segment const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); - linSys_.resWell_[seg][SPres] -= outlet_pressure.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); - } + const int seg_upwind = upwinding_segments_[seg]; + MultisegmentWellAssemble(baseif_). + assemblePressureEq(seg, seg_upwind, outlet_segment_index, + pressure_equation, outlet_pressure, linSys_); if (this->accelerationalPressureLossConsidered()) { handleAccelerationPressureLoss(seg, well_state); @@ -1535,8 +1404,8 @@ assembleICDPressureEq(const int seg, if (const auto& segment = this->segmentSet()[seg]; (segment.segmentType() == Segment::SegmentType::VALVE) && (segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve - linSys_.resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][WQTotal].value(); - linSys_.duneD_[seg][seg][SPres][WQTotal] = 1.; + MultisegmentWellAssemble(baseif_). + assembleTrivialEq(seg, this->primary_variables_evaluation_[seg][WQTotal].value(), linSys_); auto& ws = well_state.well(baseif_.indexOfWell()); ws.segments.pressure_drop_friction[seg] = 0.; @@ -1569,25 +1438,17 @@ assembleICDPressureEq(const int seg, auto& ws = well_state.well(baseif_.indexOfWell()); ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value(); - const int seg_upwind = upwinding_segments_[seg]; - linSys_.resWell_[seg][SPres] = pressure_equation.value(); - linSys_.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); - linSys_.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - linSys_.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - linSys_.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); - } - // contribution from the outlet segment const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); - linSys_.resWell_[seg][SPres] -= outlet_pressure.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); - } + const int seg_upwind = upwinding_segments_[seg]; + MultisegmentWellAssemble(baseif_). + assemblePressureEq(seg, seg_upwind, outlet_segment_index, + pressure_equation, outlet_pressure, + linSys_, + FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx), + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)); } template @@ -1623,10 +1484,10 @@ getFiniteWellResiduals(const std::vector& B_avg, for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { double residual = 0.; if (eq_idx < baseif_.numComponents()) { - residual = std::abs(linSys_.resWell_[seg][eq_idx]) * B_avg[eq_idx]; + residual = std::abs(linSys_.residual()[seg][eq_idx]) * B_avg[eq_idx]; } else { if (seg > 0) { - residual = std::abs(linSys_.resWell_[seg][eq_idx]); + residual = std::abs(linSys_.residual()[seg][eq_idx]); } } if (std::isnan(residual) || std::isinf(residual)) { @@ -1643,7 +1504,7 @@ getFiniteWellResiduals(const std::vector& B_avg, // handling the control equation residual { - const double control_residual = std::abs(linSys_.resWell_[0][numWellEq - 1]); + const double control_residual = std::abs(linSys_.residual()[0][numWellEq - 1]); if (std::isnan(control_residual) || std::isinf(control_residual)) { deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name()); return {false, residuals}; diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index 82fc4dcec..af10baf0c 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -105,15 +105,6 @@ protected: void initMatrixAndVectors(const int num_cells); void initPrimaryVariablesEvaluation() const; - void assembleControlEq(const WellState& well_state, - const GroupState& group_state, - const Schedule& schedule, - const SummaryState& summaryState, - const Well::InjectionControls& inj_controls, - const Well::ProductionControls& prod_controls, - const double rho, - DeferredLogger& deferred_logger); - void assembleDefaultPressureEq(const int seg, WellState& well_state); diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 2608313df..dd08be69b 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -19,6 +19,7 @@ */ +#include #include #include #include @@ -1562,30 +1563,18 @@ namespace Opm for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx) - segment_fluid_initial_[seg][comp_idx]) / dt; - - this->linSys_.resWell_[seg][comp_idx] += accumulation_term.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - this->linSys_.duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq); - } + MultisegmentWellAssemble(*this). + assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_); } } // considering the contributions due to flowing out from the segment { + const int seg_upwind = this->upwinding_segments_[seg]; for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { - const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_; - - const int seg_upwind = this->upwinding_segments_[seg]; - // segment_rate contains the derivatives with respect to WQTotal in seg, - // and WFrac and GFrac in seg_upwind - this->linSys_.resWell_[seg][comp_idx] -= segment_rate.value(); - this->linSys_.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - this->linSys_.duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - this->linSys_.duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); - } - // pressure derivative should be zero + const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * + this->well_efficiency_factor_; + MultisegmentWellAssemble(*this). + assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_); } } @@ -1594,19 +1583,9 @@ namespace Opm for (const int inlet : this->segment_inlets_[seg]) { for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_; - const int inlet_upwind = this->upwinding_segments_[inlet]; - // inlet_rate contains the derivatives with respect to WQTotal in inlet, - // and WFrac and GFrac in inlet_upwind - this->linSys_.resWell_[seg][comp_idx] += inlet_rate.value(); - this->linSys_.duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - this->linSys_.duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - this->linSys_.duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); - } - // pressure derivative should be zero + MultisegmentWellAssemble(*this). + assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_); } } } @@ -1647,23 +1626,8 @@ namespace Opm this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); - // subtract sum of phase fluxes in the well equations. - this->linSys_.resWell_[seg][comp_idx] += cq_s_effective.value(); - - // assemble the jacobians - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - - // also need to consider the efficiency factor when manipulating the jacobians. - this->linSys_.duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix - - // the index name for the D should be eq_idx / pv_idx - this->linSys_.duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq); - } - - for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) { - // also need to consider the efficiency factor when manipulating the jacobians. - this->linSys_.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); - } + MultisegmentWellAssemble(*this). + assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_); } } @@ -1671,13 +1635,19 @@ namespace Opm if (seg == 0) { // top segment, pressure equation is the control equation const auto& summaryState = ebosSimulator.vanguard().summaryState(); const Schedule& schedule = ebosSimulator.vanguard().schedule(); - this->assembleControlEq(well_state, + std::function gQ = [this](int a) { return this->getQs(a); }; + MultisegmentWellAssemble(*this). + assembleControlEq(well_state, group_state, schedule, summaryState, inj_controls, prod_controls, getRefDensity(), + this->getWQTotal(), + this->getBhp(), + gQ, + this->linSys_, deferred_logger); } else { const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();