Merge pull request #4306 from akva2/msw_assemble

add: MultisegmentWellAssemble

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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
+
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
+#include <opm/models/blackoil/blackoilindices.hh>
+#include <opm/models/blackoil/blackoilonephaseindices.hh>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
+#include <opm/simulators/wells/WellAssemble.hpp>
+#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
+#include <opm/simulators/wells/WellInterfaceIndices.hpp>
+
+namespace Opm {
+
+//! \brief Class administering assembler access to equation system.
+template<class Scalar, int numWellEq, int numEq>
+class MultisegmentWellEquationAccess {
+public:
+    //! \brief Constructor initializes reference to the equation system.
+    MultisegmentWellEquationAccess(MultisegmentWellEquations<Scalar,numWellEq,numEq>& eqns)
+        : eqns_(eqns)
+    {}
+
+    using BVectorWell = typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::BVectorWell;
+    using DiagMatWell = typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::DiagMatWell;
+    using OffDiatMatWell = typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::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<Scalar,numWellEq,numEq>& eqns_; //!< Reference to equation system
+};
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+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<EvalWell(const int)>& 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<EvalWell> 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<EvalWell()> 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<EvalWell()> 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<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assemblePressureLoss(const int seg,
+                     const int seg_upwind,
+                     const EvalWell& accelerationPressureLoss,
+                     Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+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<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleTrivialEq(const int seg,
+                  const Scalar value,
+                  Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
+    eqns.residual()[seg][SPres] = value;
+    eqns.D()[seg][seg][SPres][WQTotal] = 1.;
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleAccumulationTerm(const int seg,
+                         const int comp_idx,
+                         const EvalWell& accumulation_term,
+                         Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleOutflowTerm(const int seg,
+                    const int seg_upwind,
+                    const int comp_idx,
+                    const EvalWell& segment_rate,
+                    Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleInflowTerm(const int seg,
+                   const int inlet,
+                   const int inlet_upwind,
+                   const int comp_idx,
+                   const EvalWell& inlet_rate,
+                   Equations& eqns1) const
+ {
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assemblePerforationEq(const int seg,
+                      const int cell_idx,
+                      const int comp_idx,
+                      const EvalWell& cq_s_effective,
+                      Equations& eqns1) const
+{
+    MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> 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<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__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>)
+
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED
+#define OPM_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED
+
+#include <opm/input/eclipse/Schedule/Well/Well.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <functional>
+
+namespace Opm
+{
+
+class DeferredLogger;
+class GroupState;
+template<class Scalar, int numWellEq, int numEq> class MultisegmentWellEquations;
+class Schedule;
+class SummaryState;
+template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
+class WellState;
+
+//! \brief Class handling assemble of the equation system for MultisegmentWell.
+template<class FluidSystem, class Indices, class Scalar>
+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<Scalar,numWellEq,Indices::numEq>;
+    using EvalWell = DenseAd::Evaluation<Scalar, numWellEq+Indices::numEq>;
+    //! \brief Constructor initializes reference to well.
+    MultisegmentWellAssemble(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& 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<EvalWell(const int)>& 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<FluidSystem,Indices,Scalar>& well_; //!< Reference to well
+};
+
+}
+
+#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED

opm/simulators/wells/MultisegmentWellEquations.hpp

@@ -36,6 +36,7 @@ template<class M> class UMFPack;
 namespace Opm
 {
 
+template<class Scalar, int numWellEq, int numEq> class MultisegmentWellEquationAccess;
 template<class Scalar> 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<Scalar,numWellEq,numEq>;
     // 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<Dune::UMFPack<DiagMatWell> > duneDSolver_;
+    mutable std::shared_ptr<Dune::UMFPack<DiagMatWell>> duneDSolver_;
 
     // residuals of the well equations
     BVectorWell resWell_;
 
-private:
     const MultisegmentWellGeneric<Scalar>& well_; //!< Reference to well
 };
 

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -33,6 +33,7 @@
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/MSWellHelpers.hpp>
+#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/WellAssemble.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
@@ -112,7 +113,7 @@ getWellConvergence(const WellState& well_state,
     std::vector<std::vector<double>> abs_residual(this->numberOfSegments(), std::vector<double>(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<typename FluidSystem, typename Indices, typename Scalar>
-void
-MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-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<EvalWell> 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<EvalWell()> 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<EvalWell()> 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<typename FluidSystem, typename Indices, typename Scalar>
 void
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@@ -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();
+    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-    linSys_.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
+        assemblePressureLoss(seg, seg_upwind, accelerationPressureLoss, linSys_);
-    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);
-    }
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -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();
+    const int seg_upwind = upwinding_segments_[seg];
-    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-        linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
+        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();
+        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-        linSys_.duneD_[seg][seg][SPres][WQTotal] = 1.;
+            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();
+    const int seg_upwind = upwinding_segments_[seg];
-    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-        linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
+        assemblePressureEq(seg, seg_upwind, outlet_segment_index,
-    }
+                           pressure_equation, outlet_pressure,
+                           linSys_,
+                           FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
+                           FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1623,10 +1484,10 @@ getFiniteWellResiduals(const std::vector<Scalar>& 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<Scalar>& 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};

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);
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -19,6 +19,7 @@
 */
 
 
+#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
@@ -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;
-
+                    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
-                    this->linSys_.resWell_[seg][comp_idx] += accumulation_term.value();
+                        assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_);
-                    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);
-                    }
                 }
             }
             // considering the contributions due to flowing out from the segment
             {
-                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,
+                for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
-                    // and WFrac and GFrac in seg_upwind
+                    const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) *
-                    this->linSys_.resWell_[seg][comp_idx] -= segment_rate.value();
+                                                  this->well_efficiency_factor_;
-                    this->linSys_.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
+                    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
-                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                        assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_);
-                        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
                 }
             }
 
@@ -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,
+                        MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
-                        // and WFrac and GFrac in inlet_upwind
+                            assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_);
-                        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
                     }
                 }
             }
@@ -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.
+                    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
-                    this->linSys_.resWell_[seg][comp_idx] += cq_s_effective.value();
+                        assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
-
-                    // 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);
-                    }
                 }
             }
 
@@ -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<EvalWell(const int)> gQ = [this](int a) { return this->getQs(a); };
+                MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*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();