From 956916e8a05705ff0dfdc8069023e483d6c985d0 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 16 Oct 2020 15:06:14 +0200 Subject: [PATCH 1/3] WIP Add flow_dyn Adds a dynamic flow version that only takes number of equations as paramter The goal is to reduce the number of flow versions i.e reduce compile time The drawback is a slight overhead in the runtime. Initial testing indicates in the range of 2-3%. This may change.... --- CMakeLists.txt | 24 +++++ ebos/eclequilinitializer.hh | 8 +- ebos/eclfluxmodule.hh | 6 +- ebos/ecloutputblackoilmodule.hh | 3 +- ebos/eclproblem.hh | 55 ++++++----- ebos/eclwriter.hh | 4 +- flow/flow_1.cpp | 92 ++++++++++++++++++ flow/flow_1.hpp | 33 +++++++ flow/flow_2.cpp | 92 ++++++++++++++++++ flow/flow_2.hpp | 33 +++++++ flow/flow_3.cpp | 92 ++++++++++++++++++ flow/flow_3.hpp | 33 +++++++ flow/flow_4.cpp | 93 +++++++++++++++++++ flow/flow_4.hpp | 33 +++++++ opm/simulators/aquifers/AquiferInterface.hpp | 17 +--- opm/simulators/flow/BlackoilModelEbos.hpp | 16 ++-- opm/simulators/flow/Main.hpp | 42 ++++++++- .../wells/BlackoilWellModel_impl.hpp | 9 +- opm/simulators/wells/MultisegmentWell.hpp | 4 +- .../wells/MultisegmentWell_impl.hpp | 2 +- opm/simulators/wells/StandardWell.hpp | 18 ++-- opm/simulators/wells/StandardWell_impl.hpp | 46 ++++----- opm/simulators/wells/WellInterface.hpp | 17 ++-- opm/simulators/wells/WellInterface_impl.hpp | 4 +- 24 files changed, 665 insertions(+), 111 deletions(-) create mode 100644 flow/flow_1.cpp create mode 100644 flow/flow_1.hpp create mode 100644 flow/flow_2.cpp create mode 100644 flow/flow_2.hpp create mode 100644 flow/flow_3.cpp create mode 100644 flow/flow_3.hpp create mode 100644 flow/flow_4.cpp create mode 100644 flow/flow_4.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 89970ff36c9..5b67cdbff1b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,7 @@ option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON) set( USE_OPENMP_DEFAULT OFF ) # Use of OpenMP is considered experimental option(BUILD_FLOW "Build the production oriented flow simulator?" ON) option(BUILD_FLOW_BLACKOIL_ONLY "Build the production oriented flow simulator only supporting the blackoil model?" OFF) +option(BUILD_FLOW_DYN "Build the flow simulator with dynamic indicies?" OFF) option(BUILD_FLOW_VARIANTS "Build the variants for flow by default?" OFF) option(BUILD_EBOS "Build the research oriented ebos simulator?" ON) option(BUILD_EBOS_EXTENSIONS "Build the variants for various extensions of ebos by default?" OFF) @@ -293,6 +294,29 @@ opm_add_test(flow_blackoil_dunecpr DEPENDS opmsimulators LIBRARIES opmsimulators) +if (NOT BUILD_FLOW_DYN) + set(FLOW_DYN_DEFAULT_ENABLE_IF "FALSE") +else() + set(FLOW_DYN_DEFAULT_ENABLE_IF "TRUE") +endif() + +# the production oriented general-purpose ECL simulator +opm_add_test(flow_dyn + ONLY_COMPILE + ALWAYS_ENABLE + DEFAULT_ENABLE_IF ${FLOW_DYN_DEFAULT_ENABLE_IF} + DEPENDS opmsimulators + LIBRARIES opmsimulators + SOURCES + flow/flow.cpp + flow/flow_1.cpp + flow/flow_2.cpp + flow/flow_3.cpp + flow/flow_4.cpp + flow/flow_ebos_blackoil.cpp + $) +target_compile_definitions(flow_dyn PRIVATE "FLOW_DYN") + opm_add_test(flow_onephase ONLY_COMPILE DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index bfae0a64996..00aa05b68fc 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -82,7 +82,7 @@ public: FluidSystem, enableTemperature, enableEnergy, - Indices::gasEnabled, + true, enableBrine, Indices::numPhases > ScalarFluidState; @@ -116,18 +116,16 @@ public: for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (FluidSystem::phaseIsActive(phaseIdx)) fluidState.setSaturation(phaseIdx, initialState.saturation()[phaseIdx][elemIdx]); - else if (Indices::numPhases == 3) - fluidState.setSaturation(phaseIdx, 0.0); } if (FluidSystem::enableDissolvedGas()) fluidState.setRs(initialState.rs()[elemIdx]); - else if (Indices::gasEnabled) + else if (Indices::gasIsActive() && Indices::oilIsActive()) fluidState.setRs(0.0); if (FluidSystem::enableVaporizedOil()) fluidState.setRv(initialState.rv()[elemIdx]); - else if (Indices::gasEnabled) + else if (Indices::gasIsActive() && Indices::oilIsActive()) fluidState.setRv(0.0); diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index e081b682149..9995366f0da 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -106,11 +106,11 @@ class EclTransExtensiveQuantities using Evaluation = GetPropType; using GridView = GetPropType; using MaterialLaw = GetPropType; + using Indices = GetPropType; enum { dimWorld = GridView::dimensionworld }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { numPhases = FluidSystem::numPhases }; - enum { enableSolvent = getPropValue() }; enum { enableEnergy = getPropValue() }; typedef Opm::MathToolbox Toolbox; @@ -438,7 +438,7 @@ protected: volumeFlux_[phaseIdx] = pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); - if (enableSolvent && phaseIdx == gasPhaseIdx) + if (Indices::solventIsActive() && phaseIdx == gasPhaseIdx) asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); } else { @@ -456,7 +456,7 @@ protected: pressureDifference_[phaseIdx]*mob*(-transModified/faceArea); // Solvent inflow is not yet supported - if (enableSolvent && phaseIdx == gasPhaseIdx) + if (Indices::solventIsActive() && phaseIdx == gasPhaseIdx) asImp_().setSolventVolumeFlux(0.0); } } diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 9bad7555a58..ceefe4fcfa7 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -94,6 +94,7 @@ class EclOutputBlackOilModule using MaterialLawParams = GetPropType; using FluidSystem = GetPropType; using GridView = GetPropType; + using Indices = GetPropType; using Element = typename GridView::template Codim<0>::Entity; using ElementIterator = typename GridView::template Codim<0>::Iterator; @@ -388,7 +389,7 @@ public: rstKeywords["RV"] = 0; } - if (getPropValue()) + if (Indices::solventIsActive()) sSol_.resize(bufferSize, 0.0); if (getPropValue()) cPolymer_.resize(bufferSize, 0.0); diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 3cfb4562d04..5c9525c305a 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -876,6 +876,8 @@ public: initFluidSystem_(); + initBlackoilIndices_(); + // deal with DRSDT unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables(); size_t numDof = this->model().numGridDof(); @@ -1775,7 +1777,7 @@ public: values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx)); values.assignNaive(initialFluidStates_[globalDofIdx]); - if (enableSolvent) + if (Indices::solventIsActive()) values[Indices::solventSaturationIdx] = solventSaturation_[globalDofIdx]; if (enablePolymer) @@ -2180,23 +2182,23 @@ private: << "conservation as a post processing step. This is currently unsupported, " << "i.e., energy conservation is always handled fully implicitly." << std::endl; - int numDeckPhases = FluidSystem::numActivePhases(); - if (numDeckPhases < Indices::numPhases) - std::cerr << "WARNING: The number of active phases specified by the deck (" - << numDeckPhases << ") is smaller than the number of compiled-in phases (" - << Indices::numPhases << "). This usually results in a significant " - << "performance degradation compared to using a specialized simulator." << std::endl; - else if (numDeckPhases < Indices::numPhases) - throw std::runtime_error("The deck enables "+std::to_string(numDeckPhases)+" phases " - "while this simulator can only handle "+ - std::to_string(Indices::numPhases)+"."); +// int numDeckPhases = FluidSystem::numActivePhases(); +// if (numDeckPhases < Indices::numPhases) +// std::cerr << "WARNING: The number of active phases specified by the deck (" +// << numDeckPhases << ") is smaller than the number of compiled-in phases (" +// << Indices::numPhases << "). This usually results in a significant " +// << "performance degradation compared to using a specialized simulator." << std::endl; +// else if (numDeckPhases < Indices::numPhases) +// throw std::runtime_error("The deck enables "+std::to_string(numDeckPhases)+" phases " +// "while this simulator can only handle "+ +// std::to_string(Indices::numPhases)+"."); // make sure that the correct phases are active - if (FluidSystem::phaseIsActive(oilPhaseIdx) && !Indices::oilEnabled) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && !Indices::oilIsActive()) throw std::runtime_error("The deck enables oil, but this simulator cannot handle it."); - if (FluidSystem::phaseIsActive(gasPhaseIdx) && !Indices::gasEnabled) + if (FluidSystem::phaseIsActive(gasPhaseIdx) && !Indices::gasIsActive()) throw std::runtime_error("The deck enables gas, but this simulator cannot handle it."); - if (FluidSystem::phaseIsActive(waterPhaseIdx) && !Indices::waterEnabled) + if (FluidSystem::phaseIsActive(waterPhaseIdx) && !Indices::waterIsActive()) throw std::runtime_error("The deck enables water, but this simulator cannot handle it."); // the opposite cases should be fine (albeit a bit slower than what's possible) } @@ -2630,6 +2632,13 @@ private: FluidSystem::initFromState(eclState, schedule); } + void initBlackoilIndices_() + { + const auto& simulator = this->simulator(); + const auto& eclState = simulator.vanguard().eclState(); + Indices::init(eclState.runspec().phases()); + } + void readInitialCondition_() { const auto& simulator = this->simulator(); @@ -2698,7 +2707,7 @@ private: size_t numElems = this->model().numGridDof(); initialFluidStates_.resize(numElems); - if (enableSolvent) + if (Indices::solventIsActive()) solventSaturation_.resize(numElems, 0.0); if (enablePolymer) @@ -2725,13 +2734,13 @@ private: // the function and discard the unchanged result. Do not index // into 'solventSaturation_' unless solvent is enabled. { - auto ssol = enableSolvent + auto ssol = Indices::solventIsActive() ? eclWriter_->eclOutputModule().getSolventSaturation(elemIdx) : Scalar(0); processRestartSaturations_(elemFluidState, ssol); - if (enableSolvent) + if (Indices::solventIsActive()) solventSaturation_[elemIdx] = ssol; } @@ -2803,7 +2812,7 @@ private: } } - if (enableSolvent) { + if (Indices::solventIsActive()) { if (solventSaturation < smallSaturationTolerance) solventSaturation = 0.0; @@ -2818,7 +2827,7 @@ private: elemFluidState.setSaturation(phaseIdx, saturation); } } - if (enableSolvent) { + if (Indices::solventIsActive()) { solventSaturation = solventSaturation / sumSaturation; } } @@ -2934,12 +2943,12 @@ private: if (FluidSystem::enableDissolvedGas()) dofFluidState.setRs(rsData[dofIdx]); - else if (Indices::gasEnabled && Indices::oilEnabled) + else if (Indices::gasIsActive() && Indices::oilIsActive()) dofFluidState.setRs(0.0); if (FluidSystem::enableVaporizedOil()) dofFluidState.setRv(rvData[dofIdx]); - else if (Indices::gasEnabled && Indices::oilEnabled) + else if (Indices::gasIsActive() && Indices::oilIsActive()) dofFluidState.setRv(0.0); ////// @@ -2966,7 +2975,7 @@ private: const auto& eclState = vanguard.eclState(); size_t numDof = this->model().numGridDof(); - if (enableSolvent) { + if (Indices::solventIsActive()) { if (eclState.fieldProps().has_double("SSOL")) solventSaturation_ = eclState.fieldProps().get_double("SSOL"); else @@ -3149,7 +3158,7 @@ private: compIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx); break; case BCComponent::SOLVENT: - if (!enableSolvent) + if (!Indices::solventIsActive()) throw std::logic_error("solvent is disabled and you're trying to add solvent to BC"); compIdx = Indices::solventSaturationIdx; diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 8029daa6006..1a2d1e26ba4 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -166,6 +166,7 @@ class EclWriter using Scalar = GetPropType; using ElementContext = GetPropType; using FluidSystem = GetPropType; + using Indices = GetPropType; using Element = typename GridView::template Codim<0>::Entity; using ElementIterator = typename GridView::template Codim<0>::Iterator; @@ -174,7 +175,6 @@ class EclWriter typedef std::vector ScalarBuffer; enum { enableEnergy = getPropValue() }; - enum { enableSolvent = getPropValue() }; public: @@ -396,7 +396,7 @@ public: {"SWAT", Opm::UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))}, {"SGAS", Opm::UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))}, {"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy}, - {"SSOLVENT" , Opm::UnitSystem::measure::identity, enableSolvent}, + {"SSOLVENT" , Opm::UnitSystem::measure::identity, Indices::solventIsActive()}, {"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()}, {"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()}, {"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()}, diff --git a/flow/flow_1.cpp b/flow/flow_1.cpp new file mode 100644 index 00000000000..cc37752d14c --- /dev/null +++ b/flow/flow_1.cpp @@ -0,0 +1,92 @@ +/* + 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 "config.h" + +#include + +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlow1Problem { + using InheritsFrom = std::tuple; +}; +} + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef Opm::BlackOilDynIndices<1,1> type; +}; +}} + +namespace Opm { +void flow1SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig) +{ + using TypeTag = Properties::TTag::EclFlow1Problem; + using Vanguard = GetPropType; + + Vanguard::setExternalSetupTime(setupTime); + Vanguard::setExternalDeck(std::move(deck)); + Vanguard::setExternalEclState(std::move(eclState)); + Vanguard::setExternalSchedule(std::move(schedule)); + Vanguard::setExternalSummaryConfig(std::move(summaryConfig)); +} + +// ----------------- Main program ----------------- +int flow1Main(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + + Opm::FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +} diff --git a/flow/flow_1.hpp b/flow/flow_1.hpp new file mode 100644 index 00000000000..f31a7dd208b --- /dev/null +++ b/flow/flow_1.hpp @@ -0,0 +1,33 @@ +/* + 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 FLOW_1_HPP +#define FLOW_1_HPP + +#include +#include +#include +#include + +namespace Opm { +void flow1SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig); +int flow1Main(int argc, char** argv, bool outputCout, bool outputFiles); +} + +#endif // FLOW_1_HPP diff --git a/flow/flow_2.cpp b/flow/flow_2.cpp new file mode 100644 index 00000000000..8148814ae59 --- /dev/null +++ b/flow/flow_2.cpp @@ -0,0 +1,92 @@ +/* + 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 "config.h" + +#include + +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlow2Problem { + using InheritsFrom = std::tuple; +}; +} + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef Opm::BlackOilDynIndices<2,2> type; +}; +}} + +namespace Opm { +void flow2SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig) +{ + using TypeTag = Properties::TTag::EclFlow2Problem; + using Vanguard = GetPropType; + + Vanguard::setExternalSetupTime(setupTime); + Vanguard::setExternalDeck(std::move(deck)); + Vanguard::setExternalEclState(std::move(eclState)); + Vanguard::setExternalSchedule(std::move(schedule)); + Vanguard::setExternalSummaryConfig(std::move(summaryConfig)); +} + +// ----------------- Main program ----------------- +int flow2Main(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + + Opm::FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +} diff --git a/flow/flow_2.hpp b/flow/flow_2.hpp new file mode 100644 index 00000000000..b50f4ed5429 --- /dev/null +++ b/flow/flow_2.hpp @@ -0,0 +1,33 @@ +/* + 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 FLOW_2_HPP +#define FLOW_2_HPP + +#include +#include +#include +#include + +namespace Opm { +void flow2SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig); +int flow2Main(int argc, char** argv, bool outputCout, bool outputFiles); +} + +#endif // FLOW_2_HPP diff --git a/flow/flow_3.cpp b/flow/flow_3.cpp new file mode 100644 index 00000000000..e87345c6057 --- /dev/null +++ b/flow/flow_3.cpp @@ -0,0 +1,92 @@ +/* + 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 "config.h" + +#include + +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlow3Problem { + using InheritsFrom = std::tuple; +}; +} + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef Opm::BlackOilDynIndices<3,3> type; +}; +}} + +namespace Opm { +void flow3SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig) +{ + using TypeTag = Properties::TTag::EclFlow3Problem; + using Vanguard = GetPropType; + + Vanguard::setExternalSetupTime(setupTime); + Vanguard::setExternalDeck(std::move(deck)); + Vanguard::setExternalEclState(std::move(eclState)); + Vanguard::setExternalSchedule(std::move(schedule)); + Vanguard::setExternalSummaryConfig(std::move(summaryConfig)); +} + +// ----------------- Main program ----------------- +int flow3Main(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + + Opm::FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +} diff --git a/flow/flow_3.hpp b/flow/flow_3.hpp new file mode 100644 index 00000000000..533d3d4d08b --- /dev/null +++ b/flow/flow_3.hpp @@ -0,0 +1,33 @@ +/* + 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 FLOW_3_HPP +#define FLOW_3_HPP + +#include +#include +#include +#include + +namespace Opm { +void flow3SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig); +int flow3Main(int argc, char** argv, bool outputCout, bool outputFiles); +} + +#endif // FLOW_2_HPP diff --git a/flow/flow_4.cpp b/flow/flow_4.cpp new file mode 100644 index 00000000000..795c120bccd --- /dev/null +++ b/flow/flow_4.cpp @@ -0,0 +1,93 @@ +/* + 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 "config.h" + +#include + +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlow4Problem { + using InheritsFrom = std::tuple; +}; +} + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: +#warning need to make a 4-3 variant for polymer++ + typedef Opm::BlackOilDynIndices<4,4> type; +}; +}} + +namespace Opm { +void flow4SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig) +{ + using TypeTag = Properties::TTag::EclFlow4Problem; + using Vanguard = GetPropType; + + Vanguard::setExternalSetupTime(setupTime); + Vanguard::setExternalDeck(std::move(deck)); + Vanguard::setExternalEclState(std::move(eclState)); + Vanguard::setExternalSchedule(std::move(schedule)); + Vanguard::setExternalSummaryConfig(std::move(summaryConfig)); +} + +// ----------------- Main program ----------------- +int flow4Main(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + + Opm::FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +} diff --git a/flow/flow_4.hpp b/flow/flow_4.hpp new file mode 100644 index 00000000000..60363c8e8b0 --- /dev/null +++ b/flow/flow_4.hpp @@ -0,0 +1,33 @@ +/* + 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 FLOW_4_HPP +#define FLOW_4_HPP + +#include +#include +#include +#include + +namespace Opm { +void flow4SetDeck(double setupTime, std::unique_ptr deck, + std::unique_ptr eclState, + std::unique_ptr schedule, + std::unique_ptr summaryConfig); +int flow4Main(int argc, char** argv, bool outputCout, bool outputFiles); +} + +#endif // FLOW_4_HPP diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 6d9a1d25150..73dc273b7d9 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -50,24 +50,17 @@ class AquiferInterface using BlackoilIndices = GetPropType; using RateVector = GetPropType; using IntensiveQuantities = GetPropType; + using Eval = GetPropType; + using Scalar = GetPropType; + enum { enableTemperature = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableBrine = getPropValue() }; static const int numEq = BlackoilIndices::numEq; - typedef double Scalar; - - typedef DenseAd::Evaluation Eval; - - typedef Opm::BlackOilFluidState - FluidState; + + typedef typename BlackOilIntensiveQuantities::FluidState FluidState; static const auto waterCompIdx = FluidSystem::waterCompIdx; static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 66d9c6509ea..ea01ca4a2cc 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -440,17 +440,17 @@ namespace Opm { const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; Scalar pressureNew; - pressureNew = priVarsNew[Indices::pressureSwitchIdx]; + pressureNew = priVarsNew[Indices::activePressureSwitchIdx()]; Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; + saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::activeWaterSaturationIdx()]; oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { - saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx]; + saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::activeCompositionSwitchIdx()]; oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx]; } @@ -461,7 +461,7 @@ namespace Opm { const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx]; Scalar pressureOld; - pressureOld = priVarsOld[Indices::pressureSwitchIdx]; + pressureOld = priVarsOld[Indices::activePressureSwitchIdx()]; Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationOld = 1.0; @@ -473,14 +473,14 @@ namespace Opm { if (FluidSystem::numActivePhases() > 1) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; + saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::activeWaterSaturationIdx()]; oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { - saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; + saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::activeCompositionSwitchIdx()]; oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; } @@ -662,7 +662,7 @@ namespace Opm { maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue ); } - if ( has_solvent_ ) { + if ( Indices::solventIsActive() ) { B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); const auto R2 = ebosResid[cell_idx][contiSolventEqIdx]; R_sum[ contiSolventEqIdx ] += R2; @@ -801,7 +801,7 @@ namespace Opm { const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); compNames[compIdx] = FluidSystem::componentName(canonicalCompIdx); } - if (has_solvent_) { + if (Indices::solventIsActive()) { compNames[solventSaturationIdx] = "Solvent"; } if (has_polymer_) { diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp index 540d4765885..ef36dc5a828 100644 --- a/opm/simulators/flow/Main.hpp +++ b/opm/simulators/flow/Main.hpp @@ -24,7 +24,14 @@ #include -# ifndef FLOW_BLACKOIL_ONLY +# ifdef FLOW_DYN +# include +# include +# include +# include +# endif + +#if !defined(FLOW_BLACKOIL_ONLY) && !defined(FLOW_DYN) # include # include # include @@ -203,8 +210,9 @@ namespace Opm // requested. if ( false ) {} -#ifndef FLOW_BLACKOIL_ONLY // Twophase cases +#ifndef FLOW_BLACKOIL_ONLY +#ifndef FLOW_DYN else if( phases.size() == 2 ) { // oil-gas if (phases.active( Opm::Phase::GAS )) { @@ -300,7 +308,35 @@ namespace Opm std::move(summaryConfig_)); return Opm::flowEbosEnergyMain(argc_, argv_, outputCout_, outputFiles_); } -#endif // FLOW_BLACKOIL_ONLY +#endif // if NOT FLOW_DYN +#ifdef FLOW_DYN + else if( phases.size() == 1 ) { + Opm::flow1SetDeck(setupTime_, std::move(deck_), + std::move(eclipseState_), + std::move(schedule_), + std::move(summaryConfig_)); + return Opm::flow1Main(argc_, argv_, outputCout_, outputFiles_); + } else if( phases.size() == 2 ) { + Opm::flow2SetDeck(setupTime_, std::move(deck_), + std::move(eclipseState_), + std::move(schedule_), + std::move(summaryConfig_)); + return Opm::flow2Main(argc_, argv_, outputCout_, outputFiles_); + } else if ( phases.size() == 3 ) { + Opm::flow3SetDeck(setupTime_, std::move(deck_), + std::move(eclipseState_), + std::move(schedule_), + std::move(summaryConfig_)); + return Opm::flow3Main(argc_, argv_, outputCout_, outputFiles_); + } else if ( phases.size() == 4 ){ + Opm::flow4SetDeck(setupTime_, std::move(deck_), + std::move(eclipseState_), + std::move(schedule_), + std::move(summaryConfig_)); + return Opm::flow4Main(argc_, argv_, outputCout_, outputFiles_); + } +#endif // FLOW_DYN +#endif // ONLY_BLACKOIL // Blackoil case else if( phases.size() == 3 ) { Opm::flowEbosBlackoilSetDeck(setupTime_, std::move(deck_), diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 6d368c680db..372bf7e7e4b 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -86,6 +86,9 @@ namespace Opm { ebosSimulator_.model().addAuxiliaryModule(this); is_cell_perforated_.resize(local_num_cells_, false); + + has_solvent_ = Indices::solventIsActive(); + } template @@ -266,10 +269,10 @@ namespace Opm { // copy of get perfpressure in Standard well // exept for value double perf_pressure = 0.0; - if (Indices::oilEnabled) { + if (Indices::oilIsActive()) { perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); } else { - if (Indices::waterEnabled) { + if (Indices::waterIsActive()) { perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); } else { perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); @@ -1568,7 +1571,7 @@ namespace Opm { return numPhases(); } int numComp = FluidSystem::numComponents; - if (has_solvent_) { + if (Indices::solventIsActive()) { numComp ++; } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index c0694e762e8..b29be9b8d9e 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -48,7 +48,7 @@ namespace Opm /// the number of reservior equations using Base::numEq; - using Base::has_solvent; + //using Base::has_solvent; using Base::has_polymer; using Base::Water; using Base::Oil; @@ -68,7 +68,7 @@ namespace Opm static constexpr int SPres = has_gas + has_water + 1; // the number of well equations TODO: it should have a more general strategy for it - static const int numWellEq = getPropValue() ? numEq : numEq + 1; + static const int numWellEq = Indices::numWellEq; using typename Base::Scalar; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 9fd5b5bcce3..9db7661d282 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -54,7 +54,7 @@ namespace Opm , segment_phase_viscosities_(numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? { // not handling solvent or polymer for now with multisegment well - if (has_solvent) { + if (Indices::solventIsActive()) { OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet"); } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 85338ddf4fe..0af126f116f 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -72,7 +72,6 @@ namespace Opm using Base::numEq; - using Base::has_solvent; using Base::has_polymer; using Base::has_foam; using Base::has_brine; @@ -83,14 +82,8 @@ namespace Opm using BrineModule = Opm::BlackOilBrineModule; // polymer concentration and temperature are already known by the well, so - // polymer and energy conservation do not need to be considered explicitly - static const int numPolymerEq = Indices::numPolymers; - static const int numEnergyEq = Indices::numEnergy; - static const int numFoamEq = Indices::numFoam; - static const int numBrineEq = Indices::numBrine; - // number of the conservation equations - static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq - numFoamEq - numBrineEq; + static const int numWellConservationEq = Indices::numWellEq; // number of the well control equations static const int numWellControlEq = 1; // number of the well equations that will always be used @@ -106,11 +99,12 @@ namespace Opm // well control equation is always the last well equation. // TODO: in the current implementation, we use the well rate as the first primary variables for injectors, // instead of G_t. - static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0); +// static const bool gasoil = numEq == 2 &gasoil? -1000: & (Indices::compositionSwitchIdx >= 0); static const int WQTotal = 0; - static const int WFrac = gasoil? -1000: 1; - static const int GFrac = gasoil? 1: 2; - static const int SFrac = !has_solvent ? -1000 : 3; + static const int WFrac = 1; + static const int GFrac = 2; + static const int SFrac = 3; +#warning TODO make a dynamic maping // the index for Bhp in primary variables and also the index of well control equation // they both will be the last one in their respective system. // TODO: we should have indices for the well equations and well primary variables separately diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c333c740bab..84af017552a 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -200,10 +200,10 @@ namespace Opm } break; case InjectorType::GAS: - if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent + if (Indices::solventIsActive() && comp_idx == contiSolventEqIdx) { // solvent inj_frac = wsolvent(); } else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) { - inj_frac = has_solvent ? 1.0 - wsolvent() : 1.0; + inj_frac = Indices::solventIsActive() ? 1.0 - wsolvent() : 1.0; } break; case InjectorType::OIL: @@ -264,7 +264,7 @@ namespace Opm return primary_variables_evaluation_[GFrac]; } - if (has_solvent && compIdx == (unsigned)contiSolventEqIdx) { + if (Indices::solventIsActive() && compIdx == (unsigned)contiSolventEqIdx) { return primary_variables_evaluation_[SFrac]; } @@ -277,7 +277,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { well_fraction -= primary_variables_evaluation_[GFrac]; } - if (has_solvent) { + if (Indices::solventIsActive()) { well_fraction -= primary_variables_evaluation_[SFrac]; } return well_fraction; @@ -328,10 +328,10 @@ namespace Opm StandardWell::getPerfCellPressure(const typename StandardWell::FluidState& fs) const { Eval pressure; - if (Indices::oilEnabled) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { pressure = fs.pressure(FluidSystem::oilPhaseIdx); } else { - if (Indices::waterEnabled) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { pressure = fs.pressure(FluidSystem::waterPhaseIdx); } else { pressure = fs.pressure(FluidSystem::gasPhaseIdx); @@ -369,7 +369,7 @@ namespace Opm const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); b_perfcells_dense[compIdx] = extendEval(fs.invB(phaseIdx)); } - if (has_solvent) { + if (Indices::solventIsActive()) { b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); } @@ -442,7 +442,7 @@ namespace Opm volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; } - if (has_solvent) { + if (Indices::solventIsActive()) { volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx]; } @@ -624,7 +624,7 @@ namespace Opm } // Store the perforation phase flux for later usage. - if (has_solvent && componentIdx == contiSolventEqIdx) { + if (Indices::solventIsActive() && componentIdx == contiSolventEqIdx) { well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); } else { well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value(); @@ -882,7 +882,7 @@ namespace Opm const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx)); } - if (has_solvent) { + if (Indices::solventIsActive()) { mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); } } else { @@ -905,7 +905,7 @@ namespace Opm } // this may not work if viscosity and relperms has been modified? - if (has_solvent) { + if (Indices::solventIsActive()) { OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger); } } @@ -973,7 +973,7 @@ namespace Opm primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; } - if (has_solvent) { + if (Indices::solventIsActive()) { const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; @@ -1054,7 +1054,7 @@ namespace Opm } double F_solvent = 0.0; - if (has_solvent) { + if (Indices::solventIsActive()) { F_solvent = primary_variables_[SFrac]; F[pu.phase_pos[Oil]] -= F_solvent; } @@ -1064,7 +1064,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]); } - if (has_solvent) { + if (Indices::solventIsActive()) { F_solvent /= (1.0 - F[pu.phase_pos[Water]]); } F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]); @@ -1077,7 +1077,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); } - if (has_solvent) { + if (Indices::solventIsActive()) { F_solvent /= (1.0 - F[pu.phase_pos[Gas]]); } F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]); @@ -1092,7 +1092,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); } - if (has_solvent) { + if (Indices::solventIsActive()) { F_solvent /= (1.0 - F[pu.phase_pos[Oil]]); } F[pu.phase_pos[Oil]] = 0.0; @@ -1104,7 +1104,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; } - if(has_solvent) { + if(Indices::solventIsActive()) { primary_variables_[SFrac] = F_solvent; } } @@ -1138,7 +1138,7 @@ namespace Opm } double F_solvent = 0.0; - if (has_solvent) { + if (Indices::solventIsActive()) { F_solvent = primary_variables_[SFrac]; F[oil_pos] -= F_solvent; } @@ -1157,7 +1157,7 @@ namespace Opm // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. // More testing is needed to make sure this is correct for well groups and THP. - if (has_solvent){ + if (Indices::solventIsActive()){ F_solvent /= scalingFactor(contiSolventEqIdx); F[pu.phase_pos[Gas]] += F_solvent; } @@ -1935,7 +1935,7 @@ namespace Opm } // We use cell values for solvent injector - if (has_solvent) { + if (Indices::solventIsActive()) { b_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); surf_dens_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventRefDensity(); } @@ -2197,7 +2197,7 @@ namespace Opm for (int comp = 0; comp < np; ++comp) { perfRates[perf * num_components_ + comp] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(comp)]; } - if(has_solvent) { + if(Indices::solventIsActive()) { perfRates[perf * num_components_ + contiSolventEqIdx] = well_state.perfRateSolvent()[first_perf_ + perf]; } } @@ -2807,7 +2807,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ; } - if (has_solvent) { + if (Indices::solventIsActive()) { primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; } } else { // total_well_rate == 0 @@ -2825,7 +2825,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (phase == InjectorType::GAS) { primary_variables_[GFrac] = 1.0 - wsolvent(); - if (has_solvent) { + if (Indices::solventIsActive()) { primary_variables_[SFrac] = wsolvent(); } } else { diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 717fd977726..c6b18c35d54 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -82,16 +82,16 @@ namespace Opm using MaterialLaw = GetPropType; using SparseMatrixAdapter = GetPropType; using RateVector = GetPropType; + using Eval = GetPropType; + using Scalar = GetPropType; static const int numEq = Indices::numEq; - typedef double Scalar; typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BlockVector BVector; - typedef DenseAd::Evaluation Eval; - static const bool has_solvent = getPropValue(); + //static const bool has_solvent = getPropValue(); static const bool has_polymer = getPropValue(); static const bool has_energy = getPropValue(); static const bool has_temperature = getPropValue(); @@ -109,14 +109,9 @@ namespace Opm // For the conversion between the surface volume rate and reservoir voidage rate using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; - static const bool compositionSwitchEnabled = Indices::gasEnabled; - using FluidState = Opm::BlackOilFluidState; + + typedef typename BlackOilIntensiveQuantities::FluidState FluidState; + /// Constructor WellInterface(const Well& well, const int time_step, const ModelParameters& param, diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 7ff34a93fc5..4a61b544932 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -99,7 +99,7 @@ namespace Opm wsolvent_ = 0.0; - if (has_solvent && well.isInjector()) { + if (Indices::solventIsActive() && well.isInjector()) { auto injectorType = well_ecl_.injectorType(); if (injectorType == InjectorType::GAS) { wsolvent_ = well_ecl_.getSolventFraction(); @@ -1229,7 +1229,7 @@ namespace Opm return 1.0; if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx) return 0.01; - if (has_solvent && phaseIdx == contiSolventEqIdx ) + if (Indices::solventIsActive() && phaseIdx == contiSolventEqIdx ) return 0.01; // we should not come this far From ec062af506367a0a14ad69f7dc21cc61a11401c9 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 26 Oct 2020 13:16:15 +0100 Subject: [PATCH 2/3] make dynamic mapping for well models --- CMakeLists_files.cmake | 1 + flow/flow_4.cpp | 4 ++ opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 9 ++- opm/simulators/wells/MultisegmentWell.hpp | 9 +-- .../wells/MultisegmentWell_impl.hpp | 68 +++++++++++++++++-- opm/simulators/wells/StandardWell.hpp | 30 ++------ opm/simulators/wells/StandardWell_impl.hpp | 63 ++++++++++++++--- opm/simulators/wells/WellInterface.hpp | 8 ++- opm/simulators/wells/WellInterface_impl.hpp | 11 ++- 10 files changed, 156 insertions(+), 49 deletions(-) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 380613f59ae..4fe66aa7859 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -238,6 +238,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/MSWellHelpers.hpp opm/simulators/wells/BlackoilWellModel.hpp opm/simulators/wells/BlackoilWellModel_impl.hpp + opm/simulators/wells/WellIndices.hpp ) list (APPEND EXAMPLE_SOURCE_FILES diff --git a/flow/flow_4.cpp b/flow/flow_4.cpp index 795c120bccd..05d7d8b40ad 100644 --- a/flow/flow_4.cpp +++ b/flow/flow_4.cpp @@ -54,6 +54,10 @@ struct Indices #warning need to make a 4-3 variant for polymer++ typedef Opm::BlackOilDynIndices<4,4> type; }; +template +struct EnableSolvent { + static constexpr bool value = true; +}; }} namespace Opm { diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 9429febabae..dd276e5484c 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -95,7 +96,6 @@ namespace Opm { typedef typename Opm::BaseAuxiliaryModule::NeighborSet NeighborSet; static const int numEq = Indices::numEq; - static const int solventSaturationIdx = Indices::solventSaturationIdx; // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 372bf7e7e4b..91f4e9b0d6e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -338,6 +338,7 @@ namespace Opm { well_state_.disableGliftOptimization(); const int reportStepIdx = ebosSimulator_.episodeIndex(); const double simulationTime = ebosSimulator_.time(); + const Opm::Phases& phases = ebosSimulator_.vanguard().eclState().runspec().phases(); int exception_thrown = 0; try { @@ -351,7 +352,7 @@ namespace Opm { // TODO: to see whether we can postpone of the intialization of the well containers to // optimize the usage of the following several member variables for (auto& well : well_container_) { - well->init(&phase_usage_, depth_, gravity_, local_num_cells_); + well->init(phases, &phase_usage_, depth_, gravity_, local_num_cells_); } // update the updated cell flag @@ -419,6 +420,8 @@ namespace Opm { void BlackoilWellModel::wellTesting(const int timeStepIdx, const double simulationTime, Opm::DeferredLogger& deferred_logger) { const auto& wtest_config = schedule().wtestConfig(timeStepIdx); + const Opm::Phases& phases = ebosSimulator_.vanguard().eclState().runspec().phases(); + if (wtest_config.size() != 0) { // there is a WTEST request // average B factors are required for the convergence checking of well equations @@ -435,7 +438,7 @@ namespace Opm { WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger); // some preparation before the well can be used - well->init(&phase_usage_, depth_, gravity_, local_num_cells_); + well->init(phases, &phase_usage_, depth_, gravity_, local_num_cells_); const Well& wellEcl = schedule().getWell(well_name, timeStepIdx); double well_efficiency_factor = wellEcl.getEfficiencyFactor(); WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx), schedule(), timeStepIdx, well_efficiency_factor); @@ -1520,7 +1523,7 @@ namespace Opm { B += 1 / fs.invB(phaseIdx).value(); } if (has_solvent_) { - auto& B = B_avg[solventSaturationIdx]; + auto& B = B_avg[Indices::activeSolventSaturationIdx()]; B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); } } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index b29be9b8d9e..a78943325b3 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -62,11 +62,6 @@ namespace Opm static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0); static constexpr bool has_water = (Indices::waterSaturationIdx >= 0); - static constexpr int GTotal = 0; - static constexpr int WFrac = has_water ? 1: -1000; - static constexpr int GFrac = has_gas ? has_water + 1 : -1000; - static constexpr int SPres = has_gas + has_water + 1; - // the number of well equations TODO: it should have a more general strategy for it static const int numWellEq = Indices::numWellEq; @@ -107,7 +102,8 @@ namespace Opm const int first_perf_index, const std::vector& perf_data); - virtual void init(const PhaseUsage* phase_usage_arg, + virtual void init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) override; @@ -216,6 +212,7 @@ namespace Opm using Base::ebosCompIdxToFlowCompIdx; using Base::getAllowCrossFlow; using Base::scalingFactor; + using Base::wellIndices; using Base::wellIsStopped_; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 9db7661d282..ba40d1c6449 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -127,12 +127,13 @@ namespace Opm template void MultisegmentWell:: - init(const PhaseUsage* phase_usage_arg, + init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) { - Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); + Base::init(phases, phase_usage_arg, depth_arg, gravity_arg, num_cells); // TODO: for StandardWell, we need to update the perf depth here using depth_arg. // for MultisegmentWell, it is much more complicated. @@ -935,6 +936,7 @@ namespace Opm double total_seg_rate = 0.0; const int seg_index = top_segment_index + seg; // the segment pressure + const int SPres = wellIndices().activePressureIdx(); primary_variables_[seg][SPres] = well_state.segPress()[seg_index]; // TODO: under what kind of circustances, the following will be wrong? // the definition of g makes the gas phase is always the last phase @@ -942,14 +944,17 @@ namespace Opm total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p]; } + const int GTotal = wellIndices().activeTotalRateIdx(); primary_variables_[seg][GTotal] = total_seg_rate; if (std::abs(total_seg_rate) > 0.) { if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); const int water_pos = pu.phase_pos[Water]; primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate; } if (has_gas) { const int gas_pos = pu.phase_pos[Gas]; + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate; } } else { // total_seg_rate == 0 @@ -958,6 +963,7 @@ namespace Opm auto phase = well.getInjectionProperties().injectorType; if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); if (phase == InjectorType::WATER) { primary_variables_[seg][WFrac] = 1.0; } else { @@ -966,6 +972,7 @@ namespace Opm } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); if (phase == InjectorType::GAS) { primary_variables_[seg][GFrac] = 1.0; } else { @@ -975,10 +982,12 @@ namespace Opm } else if (this->isProducer()) { // producers if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); primary_variables_[seg][WFrac] = 1.0 / number_of_phases_; } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[seg][GFrac] = 1.0 / number_of_phases_; } } @@ -1109,12 +1118,14 @@ namespace Opm for (int seg = 0; seg < numberOfSegments(); ++seg) { if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; @@ -1125,11 +1136,12 @@ namespace Opm // update the segment pressure { + const int SPres = wellIndices().activePressureIdx(); const int sign = dwells[seg][SPres] > 0.? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); } - + const int GTotal = wellIndices().activeTotalRateIdx(); // update the total rate // TODO: should we have a limitation of the total rate change? { primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal]; @@ -1252,20 +1264,24 @@ namespace Opm { if (has_water && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); return primary_variables_evaluation_[seg][WFrac]; } if (has_gas && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); return primary_variables_evaluation_[seg][GFrac]; } // Oil fraction EvalWell oil_fraction = 1.0; if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); oil_fraction -= primary_variables_evaluation_[seg][WFrac]; } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); oil_fraction -= primary_variables_evaluation_[seg][GFrac]; } /* if (has_solvent) { @@ -1674,6 +1690,7 @@ namespace Opm MultisegmentWell:: getSegmentPressure(const int seg) const { + const int SPres = wellIndices().activePressureIdx(); return primary_variables_evaluation_[seg][SPres]; } @@ -1699,6 +1716,7 @@ namespace Opm getSegmentRate(const int seg, const int comp_idx) const { + const int GTotal = wellIndices().activeTotalRateIdx(); return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx); } @@ -1725,9 +1743,12 @@ namespace Opm const size_t comp_idx) const { const int seg_upwind = upwinding_segments_[seg]; + // the result will contain the derivative with resepct to GTotal in segment seg, // and the derivatives with respect to WFrac GFrac in segment seg_upwind. // the derivative with respect to SPres should be zero. + const int GTotal = wellIndices().activeTotalRateIdx(); + if (seg == 0 && this->isInjector()) { const Well& well = Base::wellEcl(); auto phase = well.getInjectionProperties().injectorType; @@ -1753,7 +1774,7 @@ namespace Opm const EvalWell segment_rate = primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx); - assert(segment_rate.derivative(SPres + numEq) == 0.); + assert(segment_rate.derivative(wellIndices().activePressureIdx() + numEq) == 0.); return segment_rate; } @@ -1767,6 +1788,7 @@ namespace Opm MultisegmentWell:: getSegmentGTotal(const int seg) const { + const int GTotal = wellIndices().activeTotalRateIdx(); return primary_variables_evaluation_[seg][GTotal]; } @@ -1916,6 +1938,7 @@ namespace Opm } // using control_eq to update the matrix and residuals + const int SPres = wellIndices().activePressureIdx(); resWell_[0][SPres] = control_eq.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq); @@ -2065,14 +2088,19 @@ namespace Opm well_state.segPressDropFriction()[seg] = friction_pressure_drop.value(); } + const int GTotal = wellIndices().activeTotalRateIdx(); + const int SPres = wellIndices().activePressureIdx(); + resWell_[seg][SPres] = pressure_equation.value(); const int seg_upwind = upwinding_segments_[seg]; duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq); duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq); if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq); } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq); } @@ -2177,13 +2205,18 @@ namespace Opm well_state.segPressDropAcceleration()[seg] = accelerationPressureLoss.value(); + const int GTotal = wellIndices().activeTotalRateIdx(); + const int SPres = wellIndices().activePressureIdx(); + resWell_[seg][SPres] -= accelerationPressureLoss.value(); duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + numEq); duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq); if (has_water) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq); } if (has_gas) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq); } } @@ -2207,12 +2240,14 @@ namespace Opm if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { const int water_pos = pu.phase_pos[Water]; + const int WFrac = wellIndices().activeWaterFractionIdx(); fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[oil_pos] -= fractions[water_pos]; } if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { const int gas_pos = pu.phase_pos[Gas]; + const int GFrac = wellIndices().activeGasFractionIdx(); fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[oil_pos] -= fractions[gas_pos]; } @@ -2250,10 +2285,12 @@ namespace Opm } if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int WFrac = wellIndices().activeWaterFractionIdx(); primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]]; } if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]]; } } @@ -2307,12 +2344,14 @@ namespace Opm if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { const int water_pos = pu.phase_pos[Water]; + const int WFrac = wellIndices().activeWaterFractionIdx(); fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[oil_pos] -= fractions[water_pos]; } if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { const int gas_pos = pu.phase_pos[Gas]; + const int GFrac = wellIndices().activeGasFractionIdx(); fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[oil_pos] -= fractions[gas_pos]; } @@ -2330,6 +2369,8 @@ namespace Opm } // calculate the phase rates based on the primary variables + const int GTotal = wellIndices().activeTotalRateIdx(); + const int SPres = wellIndices().activePressureIdx(); const double g_total = primary_variables_[seg][GTotal]; const int top_segment_index = well_state.topSegmentIndex(index_of_well_); for (int p = 0; p < number_of_phases_; ++p) { @@ -2561,12 +2602,15 @@ namespace Opm const int seg_upwind = upwinding_segments_[seg]; // segment_rate contains the derivatives with respect to GTotal in seg, // and WFrac and GFrac in seg_upwind + const int GTotal = wellIndices().activeTotalRateIdx(); resWell_[seg][comp_idx] -= segment_rate.value(); duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq); } // pressure derivative should be zero @@ -2582,12 +2626,15 @@ namespace Opm const int inlet_upwind = upwinding_segments_[inlet]; // inlet_rate contains the derivatives with respect to GTotal in inlet, // and WFrac and GFrac in inlet_upwind + const int GTotal = wellIndices().activeTotalRateIdx(); resWell_[seg][comp_idx] += inlet_rate.value(); duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq); } // pressure derivative should be zero @@ -2981,6 +3028,7 @@ namespace Opm } const double pressure_tolerance = param_.tolerance_pressure_ms_wells_; + const int SPres = wellIndices().activePressureIdx(); if (residuals[SPres] > pressure_tolerance) { sum += residuals[SPres] / pressure_tolerance; ++count; @@ -3134,6 +3182,7 @@ namespace Opm } } + const int SPres = wellIndices().activePressureIdx(); const double well_control_residual = std::abs(resWell_[0][SPres]); const int dummy_component = -1; const double max_residual_allowed = param_.max_residual_allowed_; @@ -3156,6 +3205,7 @@ namespace Opm MultisegmentWell:: updateUpwindingSegments() { + const int GTotal = wellIndices().activeTotalRateIdx(); for (int seg = 0; seg < numberOfSegments(); ++seg) { // special treatment is needed for segment 0 if (seg == 0) { @@ -3207,13 +3257,18 @@ namespace Opm well_state.segPressDropFriction()[seg] = sicd_pressure_drop.value(); const int seg_upwind = upwinding_segments_[seg]; + const int GTotal = wellIndices().activeTotalRateIdx(); + const int SPres = wellIndices().activePressureIdx(); + resWell_[seg][SPres] = pressure_equation.value(); duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq); duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq); } @@ -3252,14 +3307,17 @@ namespace Opm const auto valve_pressure_drop = pressureDropValve(seg); pressure_equation = pressure_equation - valve_pressure_drop; well_state.segPressDropFriction()[seg] = valve_pressure_drop.value(); - + const int SPres = wellIndices().activePressureIdx(); + const int GTotal = wellIndices().activeTotalRateIdx(); resWell_[seg][SPres] = pressure_equation.value(); duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq); duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq); } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 0af126f116f..1878186fdd7 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -91,24 +91,6 @@ namespace Opm static const int numStaticWellEq = numWellConservationEq + numWellControlEq; // the positions of the primary variables for StandardWell - // the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g, - // which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP. - // correspondingly, we have four well equations for blackoil model, the first three are mass - // converstation equations, and the last one is the well control equation. - // primary variables related to other components, will be before the Bhp and after F_g. - // well control equation is always the last well equation. - // TODO: in the current implementation, we use the well rate as the first primary variables for injectors, - // instead of G_t. -// static const bool gasoil = numEq == 2 &gasoil? -1000: & (Indices::compositionSwitchIdx >= 0); - static const int WQTotal = 0; - static const int WFrac = 1; - static const int GFrac = 2; - static const int SFrac = 3; -#warning TODO make a dynamic maping - // the index for Bhp in primary variables and also the index of well control equation - // they both will be the last one in their respective system. - // TODO: we should have indices for the well equations and well primary variables separately - static const int Bhp = numStaticWellEq - numWellControlEq; using typename Base::Scalar; @@ -155,7 +137,8 @@ namespace Opm const int first_perf_index, const std::vector& perf_data); - virtual void init(const PhaseUsage* phase_usage_arg, + virtual void init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) override; @@ -303,6 +286,7 @@ namespace Opm using Base::scalingFactor; using Base::scaleProductivityIndex; using Base::mostStrictBhpFromBhpLimits; + using Base::wellIndices; // protected member variables from the Base class using Base::current_step_; @@ -536,12 +520,12 @@ namespace Opm // calculate a relaxation factor to avoid overshoot of the fractions for producers // which might result in negative rates - static double relaxationFactorFractionsProducer(const std::vector& primary_variables, - const BVectorWell& dwells); + double relaxationFactorFractionsProducer(const std::vector& primary_variables, + const BVectorWell& dwells) const; // calculate a relaxation factor to avoid overshoot of total rates - static double relaxationFactorRate(const std::vector& primary_variables, - const BVectorWell& dwells); + double relaxationFactorRate(const std::vector& primary_variables, + const BVectorWell& dwells) const; virtual void wellTestingPhysical(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 84af017552a..95358dd1142 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -59,12 +59,13 @@ namespace Opm template void StandardWell:: - init(const PhaseUsage* phase_usage_arg, + init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells) { - Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); + Base::init(phases, phase_usage_arg, depth_arg, gravity_arg, num_cells); perf_depth_.resize(number_of_perforations_, 0.); for (int perf = 0; perf < number_of_perforations_; ++perf) { @@ -164,7 +165,7 @@ namespace Opm StandardWell:: getBhp() const { - return primary_variables_evaluation_[Bhp]; + return primary_variables_evaluation_[wellIndices().activePressureIdx()]; } @@ -176,6 +177,7 @@ namespace Opm StandardWell:: getWQTotal() const { + const int WQTotal = wellIndices().activeTotalRateIdx(); return primary_variables_evaluation_[WQTotal]; } @@ -190,7 +192,7 @@ namespace Opm { // Note: currently, the WQTotal definition is still depends on Injector/Producer. assert(comp_idx < num_components_); - + const int WQTotal = wellIndices().activeTotalRateIdx(); if (this->isInjector()) { // only single phase injection double inj_frac = 0.0; switch (this->wellEcl().injectorType()) { @@ -257,27 +259,33 @@ namespace Opm } if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); return primary_variables_evaluation_[WFrac]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); return primary_variables_evaluation_[GFrac]; } if (Indices::solventIsActive() && compIdx == (unsigned)contiSolventEqIdx) { + const int SFrac = wellIndices().activeSolventFractionIdx(); return primary_variables_evaluation_[SFrac]; } // Oil fraction EvalWell well_fraction(numWellEq_ + numEq, 1.0); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); well_fraction -= primary_variables_evaluation_[WFrac]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); well_fraction -= primary_variables_evaluation_[GFrac]; } if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); well_fraction -= primary_variables_evaluation_[SFrac]; } return well_fraction; @@ -378,6 +386,7 @@ namespace Opm EvalWell drawdown = pressure - well_pressure; if (this->has_polymermw && this->isInjector()) { + const int Bhp = wellIndices().activePressureIdx(); const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index]; drawdown += skin_pressure; @@ -846,6 +855,7 @@ namespace Opm // using control_eq to update the matrix and residuals // TODO: we should use a different index system for the well equations + const int Bhp = wellIndices().activePressureIdx(); resWell_[0][Bhp] = control_eq.value(); for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) { invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + numEq); @@ -958,9 +968,9 @@ namespace Opm const double relaxation_factor_fractions = (this->isProducer()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells) : 1.0; - // update the second and third well variable (The flux fractions) if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; @@ -968,12 +978,14 @@ namespace Opm } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; } if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; @@ -982,11 +994,13 @@ namespace Opm processFractions(); // updating the total rates Q_t + const int WQTotal = wellIndices().activeTotalRateIdx(); const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells); primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate; // updating the bottom hole pressure { + const int Bhp = wellIndices().activePressureIdx(); const double dBHPLimit = param_.dbhp_max_rel_; const int sign1 = dwells[0][Bhp] > 0 ? 1: -1; const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit); @@ -1014,6 +1028,7 @@ namespace Opm updateExtraPrimaryVariables(const BVectorWell& dwells) const { // for the water velocity and skin pressure + const int Bhp = wellIndices().activePressureIdx(); if (this->has_polymermw && this->isInjector()) { for (int perf = 0; perf < number_of_perforations_; ++perf) { const int wat_vel_index = Bhp + 1 + perf; @@ -1044,17 +1059,20 @@ namespace Opm F[pu.phase_pos[Oil]] = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); F[pu.phase_pos[Water]] = primary_variables_[WFrac]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]]; } double F_solvent = 0.0; if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); F_solvent = primary_variables_[SFrac]; F[pu.phase_pos[Oil]] -= F_solvent; } @@ -1099,12 +1117,15 @@ namespace Opm } if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); primary_variables_[WFrac] = F[pu.phase_pos[Water]]; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; } if(Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); primary_variables_[SFrac] = F_solvent; } } @@ -1126,6 +1147,7 @@ namespace Opm F[oil_pos] = 1.0; if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int WFrac = wellIndices().activeWaterFractionIdx(); const int water_pos = pu.phase_pos[Water]; F[water_pos] = primary_variables_[WFrac]; F[oil_pos] -= F[water_pos]; @@ -1133,12 +1155,14 @@ namespace Opm if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { const int gas_pos = pu.phase_pos[Gas]; + const int GFrac = wellIndices().activeGasFractionIdx(); F[gas_pos] = primary_variables_[GFrac]; F[oil_pos] -= F[gas_pos]; } double F_solvent = 0.0; if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); F_solvent = primary_variables_[SFrac]; F[oil_pos] -= F_solvent; } @@ -1162,10 +1186,12 @@ namespace Opm F[pu.phase_pos[Gas]] += F_solvent; } + const int Bhp = wellIndices().activePressureIdx(); well_state.bhp()[index_of_well_] = primary_variables_[Bhp]; // calculate the phase rates based on the primary variables // for producers, this is not a problem, while not sure for injectors here + const int WQTotal = wellIndices().activeTotalRateIdx(); if (this->isProducer()) { const double g_total = primary_variables_[WQTotal]; for (int p = 0; p < number_of_phases_; ++p) { @@ -2777,6 +2803,7 @@ namespace Opm // Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate // under surface condition is used here + const int WQTotal = wellIndices().activeTotalRateIdx(); if (this->isInjector()) { switch (this->wellEcl().injectorType()) { case InjectorType::WATER: @@ -2802,12 +2829,15 @@ namespace Opm if (std::abs(total_well_rate) > 0.) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ; } if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; } } else { // total_well_rate == 0 @@ -2815,6 +2845,7 @@ namespace Opm auto phase = well_ecl_.getInjectionProperties().injectorType; // only single phase injection handled if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); if (phase == InjectorType::WATER) { primary_variables_[WFrac] = 1.0; } else { @@ -2823,9 +2854,11 @@ namespace Opm } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); if (phase == InjectorType::GAS) { primary_variables_[GFrac] = 1.0 - wsolvent(); if (Indices::solventIsActive()) { + const int SFrac = wellIndices().activeSolventFractionIdx(); primary_variables_[SFrac] = wsolvent(); } } else { @@ -2839,9 +2872,11 @@ namespace Opm } else if (this->isProducer()) { // producers // TODO: the following are not addressed for the solvent case yet if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); primary_variables_[WFrac] = 1.0 / np; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); primary_variables_[GFrac] = 1.0 / np; } } else { @@ -2851,6 +2886,7 @@ namespace Opm // BHP + const int Bhp = wellIndices().activePressureIdx(); primary_variables_[Bhp] = well_state.bhp()[index_of_well_]; // other primary variables related to polymer injection @@ -3096,7 +3132,7 @@ namespace Opm double StandardWell:: relaxationFactorFractionsProducer(const std::vector& primary_variables, - const BVectorWell& dwells) + const BVectorWell& dwells) const { // TODO: not considering solvent yet // 0.95 is a experimental value, which remains to be optimized @@ -3104,11 +3140,13 @@ namespace Opm if (FluidSystem::numActivePhases() > 1) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int WFrac = wellIndices().activeWaterFractionIdx(); const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int GFrac = wellIndices().activeGasFractionIdx(); const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); } @@ -3116,6 +3154,8 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will // not be negative oil fraction later + const int GFrac = wellIndices().activeGasFractionIdx(); + const int WFrac = wellIndices().activeWaterFractionIdx(); const double original_sum = primary_variables[WFrac] + primary_variables[GFrac]; const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor; const double possible_updated_sum = original_sum - relaxed_update; @@ -3141,10 +3181,10 @@ namespace Opm double StandardWell:: relaxationFactorRate(const std::vector& primary_variables, - const BVectorWell& dwells) + const BVectorWell& dwells) const { double relaxation_factor = 1.0; - + const int WQTotal = wellIndices().activeTotalRateIdx(); // For injector, we only check the total rates to avoid sign change of rates const double original_total_rate = primary_variables[WQTotal]; const double newton_update = dwells[0][WQTotal]; @@ -3328,6 +3368,7 @@ namespace Opm updateWaterThroughput(const double dt, WellState &well_state) const { if (this->has_polymermw && this->isInjector()) { + const int Bhp = wellIndices().activePressureIdx(); for (int perf = 0; perf < number_of_perforations_; ++perf) { const double perf_water_vel = primary_variables_[Bhp + 1 + perf]; // we do not consider the formation damage due to water flowing from reservoir into wellbore @@ -3351,6 +3392,7 @@ namespace Opm std::vector& cq_s, Opm::DeferredLogger& deferred_logger) { + const int Bhp = wellIndices().activePressureIdx(); const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); const EvalWell& water_flux_s = cq_s[water_comp_idx]; const auto& fs = int_quants.fluidState(); @@ -3465,7 +3507,7 @@ namespace Opm OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); } } - + const int Bhp = wellIndices().activePressureIdx(); const double well_control_residual = std::abs(resWell_[0][Bhp]); const int dummy_component = -1; const double max_residual_allowed = param_.max_residual_allowed_; @@ -3499,6 +3541,7 @@ namespace Opm using CR = ConvergenceReport; const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance; for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int Bhp = wellIndices().activePressureIdx(); const double wat_vel_residual = res[Bhp + 1 + perf]; if (std::isnan(wat_vel_residual)) { report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, name()}); @@ -3513,6 +3556,7 @@ namespace Opm const double pskin_tol = 1000.; // 1000 pascal const auto pskin_failure_type = CR::WellFailure::Type::Pressure; for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int Bhp = wellIndices().activePressureIdx(); const double pskin_residual = res[Bhp + 1 + perf + number_of_perforations_]; if (std::isnan(pskin_residual)) { report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, name()}); @@ -3541,6 +3585,7 @@ namespace Opm // the source term related to transport of molecular weight EvalWell cq_s_polymw = cq_s_poly; if (this->isInjector()) { + const int Bhp = wellIndices().activePressureIdx(); const int wat_vel_index = Bhp + 1 + perf; const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index]; if (water_velocity > 0.) { // injecting diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index c6b18c35d54..b69c69b4484 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -145,7 +146,8 @@ namespace Opm void setGuideRate(const GuideRate* guide_rate_arg); - virtual void init(const PhaseUsage* phase_usage_arg, + virtual void init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, const int num_cells); @@ -286,6 +288,8 @@ namespace Opm protected: + WellIndices well_indices_; + // to indicate a invalid completion static const int INVALIDCOMPLETION = INT_MAX; @@ -379,6 +383,8 @@ namespace Opm double wsolvent_; + const WellIndices& wellIndices() const; + const PhaseUsage& phaseUsage() const; int flowPhaseToEbosCompIdx( const int phaseIdx ) const; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 4a61b544932..a51ce7356f8 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -123,13 +123,15 @@ namespace Opm template void WellInterface:: - init(const PhaseUsage* phase_usage_arg, + init(const Phases& phases, + const PhaseUsage* phase_usage_arg, const std::vector& /* depth_arg */, const double gravity_arg, const int /* num_cells */) { phase_usage_ = phase_usage_arg; gravity_ = gravity_arg; + well_indices_.init(phases); } @@ -254,6 +256,13 @@ namespace Opm } + template + const WellIndices& + WellInterface:: + wellIndices() const + { + return well_indices_; + } template const PhaseUsage& From ebccaa6f4d9257aa02aeefb874fe3e305f344c68 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 26 Oct 2020 14:57:03 +0100 Subject: [PATCH 3/3] add back some code --- ebos/eclproblem.hh | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 5c9525c305a..8fa71675ad1 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2182,16 +2182,16 @@ private: << "conservation as a post processing step. This is currently unsupported, " << "i.e., energy conservation is always handled fully implicitly." << std::endl; -// int numDeckPhases = FluidSystem::numActivePhases(); -// if (numDeckPhases < Indices::numPhases) -// std::cerr << "WARNING: The number of active phases specified by the deck (" -// << numDeckPhases << ") is smaller than the number of compiled-in phases (" -// << Indices::numPhases << "). This usually results in a significant " -// << "performance degradation compared to using a specialized simulator." << std::endl; -// else if (numDeckPhases < Indices::numPhases) -// throw std::runtime_error("The deck enables "+std::to_string(numDeckPhases)+" phases " -// "while this simulator can only handle "+ -// std::to_string(Indices::numPhases)+"."); + int numDeckPhases = FluidSystem::numActivePhases(); + if (numDeckPhases < Indices::numPhases) + std::cerr << "WARNING: The number of active phases specified by the deck (" + << numDeckPhases << ") is smaller than the number of compiled-in phases (" + << Indices::numPhases << "). This usually results in a significant " + << "performance degradation compared to using a specialized simulator." << std::endl; + else if (numDeckPhases < Indices::numPhases) + throw std::runtime_error("The deck enables "+std::to_string(numDeckPhases)+" phases " + "while this simulator can only handle "+ + std::to_string(Indices::numPhases)+"."); // make sure that the correct phases are active if (FluidSystem::phaseIsActive(oilPhaseIdx) && !Indices::oilIsActive())