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