From 4ebde4e00382f30550f25b44b302a570a31cddb9 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 11:57:37 +0100 Subject: [PATCH 01/12] added: MultisegmentWellAssemble this handles assembly of the equation system for multisegment well. start by moving assembleControlEq into the new class --- CMakeLists_files.cmake | 2 + .../wells/MultisegmentWellAssemble.cpp | 186 ++++++++++++++++++ .../wells/MultisegmentWellAssemble.hpp | 76 +++++++ opm/simulators/wells/MultisegmentWellEval.cpp | 114 ----------- opm/simulators/wells/MultisegmentWellEval.hpp | 9 - .../wells/MultisegmentWell_impl.hpp | 10 +- 6 files changed, 273 insertions(+), 124 deletions(-) create mode 100644 opm/simulators/wells/MultisegmentWellAssemble.cpp create mode 100644 opm/simulators/wells/MultisegmentWellAssemble.hpp 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..114f0f6ed --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -0,0 +1,186 @@ +/* + 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 { + +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 int SPres, + const std::function& getQs, + Equations& eqns, + 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); + } + + // using control_eq to update the matrix and residuals + eqns.resWell_[0][SPres] = control_eq.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); + } +} + +#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,true,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,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..a09349e01 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -0,0 +1,76 @@ +/* + 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 +{ +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 int SPres, + const std::function& getQs, + Equations& eqns, + DeferredLogger& deferred_logger) const; + +private: + const WellInterfaceIndices& well_; //!< Reference to well +}; + +} + +#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 2514ef879..5da357eca 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1107,120 +1107,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:: 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..0418911d4 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -19,6 +19,7 @@ */ +#include #include #include #include @@ -1671,13 +1672,20 @@ 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(), + SPres, + gQ, + this->linSys_, deferred_logger); } else { const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem(); From d64508f3b8e3be4aa2eb749f1bf053040f18d34f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 02/12] added: MultisegmentWellAssemble::assemblePressureLoss extracted from MultisegmentWellEval::handleAccelerationPressureLoss --- .../wells/MultisegmentWellAssemble.cpp | 21 ++++++++++++++++- .../wells/MultisegmentWellAssemble.hpp | 23 ++++++++++++++++++- opm/simulators/wells/MultisegmentWellEval.cpp | 12 +++------- .../wells/MultisegmentWell_impl.hpp | 1 - 4 files changed, 45 insertions(+), 12 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 114f0f6ed..3b16ef468 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -48,7 +48,6 @@ assembleControlEq(const WellState& well_state, const double rho, const EvalWell& wqTotal, const EvalWell& bhp, - const int SPres, const std::function& getQs, Equations& eqns, DeferredLogger& deferred_logger) const @@ -155,6 +154,24 @@ assembleControlEq(const WellState& well_state, } } +template +void MultisegmentWellAssemble:: +assemblePressureLoss(const int seg, + const int seg_upwind, + const EvalWell& accelerationPressureLoss, + Equations& eqns) const +{ + eqns.resWell_[seg][SPres] -= accelerationPressureLoss.value(); + eqns.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); + eqns.duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq); + if constexpr (has_wfrac_variable) { + eqns.duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); + } + if constexpr (has_gfrac_variable) { + eqns.duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); + } +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; @@ -175,8 +192,10 @@ 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>) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index a09349e01..bc88d9247 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -43,6 +43,21 @@ class WellState; 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; @@ -62,11 +77,17 @@ public: const double rho, const EvalWell& wqTotal, const EvalWell& bhp, - const int SPres, 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 5da357eca..a978216de 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -1148,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 diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 0418911d4..ec492ae12 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1683,7 +1683,6 @@ namespace Opm getRefDensity(), this->getWQTotal(), this->getBhp(), - SPres, gQ, this->linSys_, deferred_logger); From 05a4ca85a79e1542cce7b08fd3da2f2e389a83e7 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 03/12] added: MultisegmentWellAssemble::assemblePressureEq extracted from MultisegmentWellEval::assembleDefaultPressureEq --- .../wells/MultisegmentWellAssemble.cpp | 28 +++++++++++++++++++ .../wells/MultisegmentWellAssemble.hpp | 10 +++++++ opm/simulators/wells/MultisegmentWellEval.cpp | 19 +++---------- 3 files changed, 42 insertions(+), 15 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 3b16ef468..78cd2be51 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -172,6 +172,34 @@ assemblePressureLoss(const int seg, } } +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& eqns, + bool wfrac, + bool gfrac) const +{ + eqns.resWell_[seg][SPres] = pressure_equation.value(); + eqns.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); + eqns.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); + if (wfrac) { + eqns.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); + } + if (gfrac) { + eqns.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); + } + + // contribution from the outlet segment + eqns.resWell_[seg][SPres] -= outlet_pressure.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); + } +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index bc88d9247..4ded47a60 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -88,6 +88,16 @@ public: 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index a978216de..c8b4643fe 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1180,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); From 32dce644d36981f91d7f2cf972c3bb6b0b386f7d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 04/12] added: MultisegmentWellAssemble::assembleTrivialEq extracted from MultisegmentWellEval::assembleICDPressureEq --- opm/simulators/wells/MultisegmentWellAssemble.cpp | 10 ++++++++++ opm/simulators/wells/MultisegmentWellAssemble.hpp | 5 +++++ opm/simulators/wells/MultisegmentWellEval.cpp | 4 ++-- 3 files changed, 17 insertions(+), 2 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 78cd2be51..a281e15f1 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -200,6 +200,16 @@ assemblePressureEq(const int seg, } } +template +void MultisegmentWellAssemble:: +assembleTrivialEq(const int seg, + const Scalar value, + Equations& eqns) const +{ + eqns.resWell_[seg][SPres] = value; + eqns.duneD_[seg][seg][SPres][WQTotal] = 1.; +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 4ded47a60..438a7b71a 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -98,6 +98,11 @@ public: 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index c8b4643fe..85ed57946 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1404,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.; From 1952ca1e5c7d98e6ddc59018cfc2aaadc266af3e Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 05/12] use MultisegmentWellAssemble::assemblePressureEq in assembleICDPressureEq --- opm/simulators/wells/MultisegmentWellEval.cpp | 22 ++++++------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 85ed57946..ebbf30854 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1438,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 From 517dd497713df95b2a13ab1ce804db7d9f012dc5 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 06/12] added: MultisegmentWellAssemble::assembleAccumulationTerm extracted from MultisegmentWell::assembleWellEqWithoutIteration --- opm/simulators/wells/MultisegmentWellAssemble.cpp | 13 +++++++++++++ opm/simulators/wells/MultisegmentWellAssemble.hpp | 6 ++++++ opm/simulators/wells/MultisegmentWell_impl.hpp | 7 ++----- 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index a281e15f1..95124d9ed 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -210,6 +210,19 @@ assembleTrivialEq(const int seg, eqns.duneD_[seg][seg][SPres][WQTotal] = 1.; } +template +void MultisegmentWellAssemble:: +assembleAccumulationTerm(const int seg, + const int comp_idx, + const EvalWell& accumulation_term, + Equations& eqns) const +{ + eqns.resWell_[seg][comp_idx] += accumulation_term.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + eqns.duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq); + } +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 438a7b71a..704f2c7dc 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -103,6 +103,12 @@ public: 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ec492ae12..56ccee4d3 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1563,11 +1563,8 @@ 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 From d5bbccde650fe785a3f54a7a649152c088ce635f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 07/12] added: MultisegmentWellAssemble::assembleOutflowTerm extracted from MultisegmentWell::assembleWellEqWithoutIteration --- .../wells/MultisegmentWellAssemble.cpp | 19 +++++++++++++++++++ .../wells/MultisegmentWellAssemble.hpp | 7 +++++++ .../wells/MultisegmentWell_impl.hpp | 19 +++++-------------- 3 files changed, 31 insertions(+), 14 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 95124d9ed..827aa1e0b 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -223,6 +223,25 @@ assembleAccumulationTerm(const int seg, } } +template +void MultisegmentWellAssemble:: +assembleOutflowTerm(const int seg, + const int seg_upwind, + const int comp_idx, + const EvalWell& segment_rate, + Equations& eqns) const +{ + eqns.resWell_[seg][comp_idx] -= segment_rate.value(); + eqns.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + eqns.duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + eqns.duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); + } + // pressure derivative should be zero +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 704f2c7dc..6eb246d44 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -109,6 +109,13 @@ public: 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 56ccee4d3..480fbcc43 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1569,21 +1569,12 @@ namespace Opm } // 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_); } } From b1d1e47e28530c19239b3db7e28f11921f773046 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 08/12] added: MultisegmentWellAssemble::assembleInflowTerm extracted from MultisegmentWell::assembleWellEqWithoutIteration --- .../wells/MultisegmentWellAssemble.cpp | 20 +++++++++++++++++++ .../wells/MultisegmentWellAssemble.hpp | 8 ++++++++ .../wells/MultisegmentWell_impl.hpp | 14 ++----------- 3 files changed, 30 insertions(+), 12 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 827aa1e0b..f6a7c7254 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -242,6 +242,26 @@ assembleOutflowTerm(const int seg, // 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& eqns) const + { + eqns.resWell_[seg][comp_idx] += inlet_rate.value(); + eqns.duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + eqns.duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + eqns.duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); + } + // pressure derivative should be zero +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 6eb246d44..43e768614 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -116,6 +116,14 @@ public: 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; + private: const WellInterfaceIndices& well_; //!< Reference to well }; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 480fbcc43..9e1e3f97a 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1583,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_); } } } From 6011a42246fc17f5314a5b6126d5b17c76095dbc Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:09:43 +0100 Subject: [PATCH 09/12] added: MultisegmentWellAssemble::assemblePerforationEq extracted from MultisegmentWell::assembleWellEqWithoutIteration --- .../wells/MultisegmentWellAssemble.cpp | 26 +++++++++++++++++++ .../wells/MultisegmentWellAssemble.hpp | 7 +++++ .../wells/MultisegmentWell_impl.hpp | 19 ++------------ 3 files changed, 35 insertions(+), 17 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index f6a7c7254..7ff23370e 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -262,6 +262,32 @@ assembleInflowTerm(const int seg, // 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& eqns) const +{ + // subtract sum of phase fluxes in the well equations. + eqns.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. + eqns.duneC_[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.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. + eqns.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); + } +} + #define INSTANCE(...) \ template class MultisegmentWellAssemble,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 43e768614..693aa0c0e 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -124,6 +124,13 @@ public: 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 }; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 9e1e3f97a..dd08be69b 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1626,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_); } } From 45457613741f0740ce54907868fac2391809d5ea Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:50:19 +0100 Subject: [PATCH 10/12] added: MultisegmentWellEquationAccess this is a proxy class for accessing the equation system in MultisegmentWellAssemble. use the new class for vector/matrix access in MultisegmentWellAssemble. --- .../wells/MultisegmentWellAssemble.cpp | 125 ++++++++++++------ 1 file changed, 87 insertions(+), 38 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index 7ff23370e..5a8382bfd 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -37,6 +37,47 @@ 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, @@ -49,7 +90,7 @@ assembleControlEq(const WellState& well_state, const EvalWell& wqTotal, const EvalWell& bhp, const std::function& getQs, - Equations& eqns, + Equations& eqns1, DeferredLogger& deferred_logger) const { static constexpr int Gas = BlackoilPhases::Vapour; @@ -147,10 +188,11 @@ assembleControlEq(const WellState& well_state, deferred_logger); } + MultisegmentWellEquationAccess eqns(eqns1); // using control_eq to update the matrix and residuals - eqns.resWell_[0][SPres] = control_eq.value(); + eqns.residual()[0][SPres] = control_eq.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - eqns.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); + eqns.D()[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq); } } @@ -159,16 +201,17 @@ void MultisegmentWellAssemble:: assemblePressureLoss(const int seg, const int seg_upwind, const EvalWell& accelerationPressureLoss, - Equations& eqns) const + Equations& eqns1) const { - eqns.resWell_[seg][SPres] -= accelerationPressureLoss.value(); - eqns.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); - eqns.duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq); + 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.duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq); } if constexpr (has_gfrac_variable) { - eqns.duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq); } } @@ -179,24 +222,25 @@ assemblePressureEq(const int seg, const int outlet_segment_index, const EvalWell& pressure_equation, const EvalWell& outlet_pressure, - Equations& eqns, + Equations& eqns1, bool wfrac, bool gfrac) const { - eqns.resWell_[seg][SPres] = pressure_equation.value(); - eqns.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); - eqns.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq); + 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.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq); } if (gfrac) { - eqns.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq); } // contribution from the outlet segment - eqns.resWell_[seg][SPres] -= outlet_pressure.value(); + eqns.residual()[seg][SPres] -= outlet_pressure.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - eqns.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); + eqns.D()[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq); } } @@ -204,10 +248,11 @@ template void MultisegmentWellAssemble:: assembleTrivialEq(const int seg, const Scalar value, - Equations& eqns) const + Equations& eqns1) const { - eqns.resWell_[seg][SPres] = value; - eqns.duneD_[seg][seg][SPres][WQTotal] = 1.; + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][SPres] = value; + eqns.D()[seg][seg][SPres][WQTotal] = 1.; } template @@ -215,11 +260,12 @@ void MultisegmentWellAssemble:: assembleAccumulationTerm(const int seg, const int comp_idx, const EvalWell& accumulation_term, - Equations& eqns) const + Equations& eqns1) const { - eqns.resWell_[seg][comp_idx] += accumulation_term.value(); + MultisegmentWellEquationAccess eqns(eqns1); + eqns.residual()[seg][comp_idx] += accumulation_term.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - eqns.duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq); + eqns.D()[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq); } } @@ -229,15 +275,16 @@ assembleOutflowTerm(const int seg, const int seg_upwind, const int comp_idx, const EvalWell& segment_rate, - Equations& eqns) const + Equations& eqns1) const { - eqns.resWell_[seg][comp_idx] -= segment_rate.value(); - eqns.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq); + 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.duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - eqns.duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); + eqns.D()[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); } // pressure derivative should be zero } @@ -249,15 +296,16 @@ assembleInflowTerm(const int seg, const int inlet_upwind, const int comp_idx, const EvalWell& inlet_rate, - Equations& eqns) const + Equations& eqns1) const { - eqns.resWell_[seg][comp_idx] += inlet_rate.value(); - eqns.duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq); + 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.duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); + eqns.D()[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - eqns.duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); + eqns.D()[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); } // pressure derivative should be zero } @@ -268,23 +316,24 @@ assemblePerforationEq(const int seg, const int cell_idx, const int comp_idx, const EvalWell& cq_s_effective, - Equations& eqns) const + Equations& eqns1) const { + MultisegmentWellEquationAccess eqns(eqns1); // subtract sum of phase fluxes in the well equations. - eqns.resWell_[seg][comp_idx] += cq_s_effective.value(); + 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.duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // input in transformed matrix + 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.duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq); + 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.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); + eqns.B()[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); } } From 2d154b50bb1d2762b77961c67b5f870eba332567 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 12:59:16 +0100 Subject: [PATCH 11/12] added: MultisegmentWellEquations::residual() this returns a const reference to the residual vector. use this in MultisegmentWellEval --- opm/simulators/wells/MultisegmentWellEquations.hpp | 6 ++++++ opm/simulators/wells/MultisegmentWellEval.cpp | 10 +++++----- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index 9bdb79d0a..34ba7e631 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -110,6 +110,12 @@ 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_; + } + // two off-diagonal matrices OffDiagMatWell duneB_; OffDiagMatWell duneC_; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index ebbf30854..3db453e69 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -113,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]); } } @@ -178,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); @@ -1484,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)) { @@ -1504,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}; From b102103e268b46d5e93bb083e4fc94a13519e01c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 18 Nov 2022 13:00:01 +0100 Subject: [PATCH 12/12] changed: make MultisegmentWellEquations data members private use a friend declaration for MultisegmentWellEquationAccess to give the assembler access to the matrices/vectors. --- opm/simulators/wells/MultisegmentWellEquations.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index 34ba7e631..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; @@ -116,6 +117,8 @@ public: return resWell_; } + private: + friend class MultisegmentWellEquationAccess; // two off-diagonal matrices OffDiagMatWell duneB_; OffDiagMatWell duneC_; @@ -125,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 };