diff --git a/CMakeLists.txt b/CMakeLists.txt index 65eea1ad0..231dc9bab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -634,6 +634,7 @@ opm_add_test(flowexp_blackoil set(FLOWEXP_COMPONENTS_SOURCES) foreach(component IN LISTS OPM_COMPILE_COMPONENTS) list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}.cpp) + list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}_2p.cpp) endforeach() opm_add_test(flowexp_comp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 51372468d..28d3fefa9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/BlackoilModelParameters.cpp opm/simulators/flow/BlackoilModelConvergenceMonitor.cpp opm/simulators/flow/CollectDataOnIORank.cpp + opm/simulators/flow/CompositionalContainer.cpp opm/simulators/flow/ConvergenceOutputConfiguration.cpp opm/simulators/flow/EclGenericWriter.cpp opm/simulators/flow/ExtboContainer.cpp @@ -119,6 +120,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/SimulatorReportBanners.cpp opm/simulators/flow/SimulatorSerializer.cpp opm/simulators/flow/SolutionContainers.cpp + opm/simulators/flow/TracerContainer.cpp opm/simulators/flow/Transmissibility.cpp opm/simulators/flow/ValidationFunctions.cpp opm/simulators/flow/equil/EquilibrationHelpers.cpp @@ -821,6 +823,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/BlackoilModelProperties.hpp opm/simulators/flow/CollectDataOnIORank.hpp opm/simulators/flow/CollectDataOnIORank_impl.hpp + opm/simulators/flow/CompositionalContainer.hpp opm/simulators/flow/ConvergenceOutputConfiguration.hpp opm/simulators/flow/countGlobalCells.hpp opm/simulators/flow/CpGridVanguard.hpp @@ -878,6 +881,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/SolutionContainers.hpp opm/simulators/flow/SubDomain.hpp opm/simulators/flow/TTagFlowProblemTPFA.hpp + opm/simulators/flow/TracerContainer.hpp opm/simulators/flow/TracerModel.hpp opm/simulators/flow/Transmissibility.hpp opm/simulators/flow/Transmissibility_impl.hpp diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index 789e2c40a..07eda2d92 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -38,7 +38,8 @@ #include #include #include -#include +#include +#include #include #include #include @@ -76,6 +77,14 @@ struct NumComp { static constexpr int value = 3; }; +template +struct EnableDummyWater { using type = UndefinedProperty; }; + +template +struct EnableDummyWater { + static constexpr bool value = true; +}; + // Set the grid type: --->2D template struct Grid { using type = Dune::YaspGrid; }; @@ -104,9 +113,10 @@ struct FluidSystem private: using Scalar = GetPropType; static constexpr int num_comp = getPropValue(); + static constexpr bool enable_water = getPropValue(); public: - using type = Opm::GenericOilGasFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; // Set the material Law @@ -118,8 +128,8 @@ private: enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; using Scalar = GetPropType; - using Traits = Opm::TwoPhaseMaterialTraits; @@ -195,6 +205,8 @@ class CO2PTProblem : public GetPropType enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { pressure0Idx = Indices::pressure0Idx }; + enum { z0Idx = Indices::z0Idx }; enum { numComponents = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; @@ -237,6 +249,21 @@ public: porosity_ = 0.1; } + void initWaterPVT() + { + using WaterPvt = typename FluidSystem::WaterPvt; + std::shared_ptr waterPvt = std::make_shared(); + waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater); + auto& ccWaterPvt = waterPvt->template getRealPvt(); + ccWaterPvt.setNumRegions(/*numPvtRegions=*/1); + Scalar rhoRefW = 1037.0; // [kg] + ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, /*rhoRefO=*/Scalar(0.0), /*rhoRefG=*/Scalar(0.0), rhoRefW); + ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4); + ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10); + waterPvt->initEnd(); + FluidSystem::setWaterPvt(waterPvt); + } + template const DimVector& gravity([[maybe_unused]]const Context& context, @@ -262,8 +289,12 @@ public: void finishInit() { ParentType::finishInit(); + // initialize fixed parameters; temperature, permeability, porosity initPetrophysics(); + + // Initialize water pvt + initWaterPVT(); } /*! @@ -358,8 +389,8 @@ public: // Calculate storage terms PrimaryVariables storageL, storageG; - this->model().globalPhaseStorage(storageL, /*phaseIdx=*/0); - this->model().globalPhaseStorage(storageG, /*phaseIdx=*/1); + this->model().globalPhaseStorage(storageL, /*phaseIdx=*/oilPhaseIdx); + this->model().globalPhaseStorage(storageG, /*phaseIdx=*/gasPhaseIdx); // Write mass balance information for rank 0 // if (this->gridView().comm().rank() == 0) { @@ -456,12 +487,12 @@ private: int prod = Parameters::Get() - 1; int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx); ComponentVector comp; - comp[0] = Evaluation::createVariable(0.5, 1); - comp[1] = Evaluation::createVariable(0.3, 2); + comp[0] = Evaluation::createVariable(0.5, z0Idx); + comp[1] = Evaluation::createVariable(0.3, z0Idx + 1); comp[2] = 1. - comp[0] - comp[1]; if (spatialIdx == inj) { - comp[0] = Evaluation::createVariable(0.99, 1); - comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2); + comp[0] = Evaluation::createVariable(0.99, z0Idx); + comp[1] = Evaluation::createVariable(0.01 - 1e-3, z0Idx + 1); comp[2] = 1. - comp[0] - comp[1]; } @@ -475,10 +506,11 @@ private: if (spatialIdx == prod) { p0 *= 0.5; } - Evaluation p_init = Evaluation::createVariable(p0, 0); + Evaluation p_init = Evaluation::createVariable(p0, pressure0Idx); fs.setPressure(FluidSystem::oilPhaseIdx, p_init); fs.setPressure(FluidSystem::gasPhaseIdx, p_init); + fs.setPressure(FluidSystem::waterPhaseIdx, p_init); fs.setTemperature(temperature_); diff --git a/flowexperimental/comp/flowexp_comp.cpp b/flowexperimental/comp/flowexp_comp.cpp index 878811398..afda32676 100644 --- a/flowexperimental/comp/flowexp_comp.cpp +++ b/flowexperimental/comp/flowexp_comp.cpp @@ -36,10 +36,13 @@ template std::tuple -runComponent(int runtimeComponent, int argc, char** argv) +runComponent(int runtimeComponent, bool water, int argc, char** argv) { if (runtimeComponent == compileTimeComponent) { - return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + if (water) + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); } return std::make_tuple(false, EXIT_FAILURE); } @@ -63,18 +66,21 @@ runComponent(int runtimeComponent, int argc, char** argv) */ template std::tuple -runComponent(int runtimecomponent, int argc, char** argv) +runComponent(int runtimecomponent, bool water, int argc, char** argv) { if (currentCompileTimeComponent == runtimecomponent) { - return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + if (water) + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); } - return runComponent(runtimecomponent, argc, argv); + return runComponent(runtimecomponent, water, argc, argv); } int main(int argc, char** argv) { - using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0>; + using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0, true>; Opm::registerEclTimeSteppingParameters(); // At the moment, this is probably as optimal as can be. @@ -90,9 +96,11 @@ main(int argc, char** argv) = Opm::Parser {}.parseFile(inputFilename, Opm::ParseContext {}, std::vector {Opm::Ecl::SectionType::RUNSPEC}); const auto runspec = Opm::Runspec(deck); const auto numComps = runspec.numComps(); + const auto& phases = runspec.phases(); + const auto wat = phases.active(Opm::Phase::WATER); auto [componentSupported, executionStatus] - = runComponent(numComps, argc, argv); + = runComponent(numComps, wat, argc, argv); if (!componentSupported) { fmt::print("Deck has {} components, not supported. In this build of the simulator, we support the " diff --git a/flowexperimental/comp/flowexp_comp.hpp b/flowexperimental/comp/flowexp_comp.hpp index 213d63c07..02920258e 100644 --- a/flowexperimental/comp/flowexp_comp.hpp +++ b/flowexperimental/comp/flowexp_comp.hpp @@ -20,7 +20,7 @@ #define FLOWEXP_COMP_HPP #include -#include +#include #include #include @@ -39,7 +39,7 @@ // suggestTimeStep is taken from newton solver in problem.limitTimestep namespace Opm { -template +template int dispatchFlowExpComp(int argc, char** argv); } @@ -47,15 +47,15 @@ int dispatchFlowExpComp(int argc, char** argv); namespace Opm::Properties { namespace TTag { -template +template struct FlowExpCompProblem { using InheritsFrom = std::tuple; }; } -template -struct SparseMatrixAdapter> +template +struct SparseMatrixAdapter> { private: using Scalar = GetPropType; @@ -109,30 +109,30 @@ struct LocalLinearizerSplice #endif // Set the problem property -template -struct Problem> +template +struct Problem> { using type = FlowProblemComp; }; -template -struct AquiferModel> { +template +struct AquiferModel> { using type = EmptyModel; }; -template -struct WellModel> { +template +struct WellModel> { using type = EmptyModel; }; -template -struct TracerModel> { +template +struct TracerModel> { using type = EmptyModel; }; -template -struct FlashSolver> { +template +struct FlashSolver> { private: using Scalar = GetPropType; using FluidSystem = GetPropType; @@ -147,10 +147,18 @@ template struct NumComp { using type = UndefinedProperty; }; // TODO: this is unfortunate, have to check why we need to hard-code it -template -struct NumComp> { +template +struct NumComp> { static constexpr int value = NumComp_; }; + +template +struct EnableDummyWater { using type = UndefinedProperty; }; + +template +struct EnableDummyWater> { + static constexpr bool value = EnableWater_; +}; #if 0 struct Temperature { using type = UndefinedProperty; }; @@ -161,26 +169,29 @@ struct Temperature { using type = UndefinedProperty; }; }; #endif -template -struct FluidSystem> +template +struct FluidSystem> { private: using Scalar = GetPropType; static constexpr int num_comp = getPropValue(); + static constexpr bool enable_water = getPropValue(); public: - using type = Opm::GenericOilGasFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; -template -struct EnableMech> { +template +struct EnableMech> { static constexpr bool value = false; }; -template -struct EnableDisgasInWater> { static constexpr bool value = false; }; +template +struct EnableDisgasInWater> { + static constexpr bool value = false; +}; -template -struct Stencil> +template +struct Stencil> { private: using Scalar = GetPropType; @@ -190,62 +201,62 @@ public: using type = EcfvStencil; }; -template -struct EnableApiTracking> { +template +struct EnableApiTracking> { static constexpr bool value = false; }; -template -struct EnableTemperature> { +template +struct EnableTemperature> { static constexpr bool value = false; }; -template -struct EnableSaltPrecipitation> { +template +struct EnableSaltPrecipitation> { static constexpr bool value = false; }; -template -struct EnablePolymerMW> { +template +struct EnablePolymerMW> { static constexpr bool value = false; }; -template -struct EnablePolymer> { +template +struct EnablePolymer> { static constexpr bool value = false; }; -template -struct EnableDispersion> { +template +struct EnableDispersion> { static constexpr bool value = false; }; -template -struct EnableBrine> { +template +struct EnableBrine> { static constexpr bool value = false; }; -template -struct EnableVapwat> { +template +struct EnableVapwat> { static constexpr bool value = false; }; -template -struct EnableSolvent> { +template +struct EnableSolvent> { static constexpr bool value = false; }; -template -struct EnableEnergy> { +template +struct EnableEnergy> { static constexpr bool value = false; }; -template -struct EnableFoam> { +template +struct EnableFoam> { static constexpr bool value = false; }; -template -struct EnableExtbo> { +template +struct EnableExtbo> { static constexpr bool value = false; }; -template -struct EnableMICP> { +template +struct EnableMICP> { static constexpr bool value = false; }; diff --git a/flowexperimental/comp/flowexp_comp2.cpp b/flowexperimental/comp/flowexp_comp2.cpp index 6886e0870..1da968da6 100644 --- a/flowexperimental/comp/flowexp_comp2.cpp +++ b/flowexperimental/comp/flowexp_comp2.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<2>(int argc, char** argv) +int dispatchFlowExpComp<2, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp2_2p.cpp b/flowexperimental/comp/flowexp_comp2_2p.cpp new file mode 100644 index 000000000..984bdb1a7 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp2_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<2, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp3.cpp b/flowexperimental/comp/flowexp_comp3.cpp index 0e5725092..0afadad14 100644 --- a/flowexperimental/comp/flowexp_comp3.cpp +++ b/flowexperimental/comp/flowexp_comp3.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<3>(int argc, char** argv) +int dispatchFlowExpComp<3, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp3_2p.cpp b/flowexperimental/comp/flowexp_comp3_2p.cpp new file mode 100644 index 000000000..d89d51533 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp3_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<3, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp4.cpp b/flowexperimental/comp/flowexp_comp4.cpp index ddb834075..7c018af27 100644 --- a/flowexperimental/comp/flowexp_comp4.cpp +++ b/flowexperimental/comp/flowexp_comp4.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<4>(int argc, char** argv) +int dispatchFlowExpComp<4, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp4_2p.cpp b/flowexperimental/comp/flowexp_comp4_2p.cpp new file mode 100644 index 000000000..68446bc24 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp4_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<4, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp5.cpp b/flowexperimental/comp/flowexp_comp5.cpp index a92539b26..d7b9cf5b7 100644 --- a/flowexperimental/comp/flowexp_comp5.cpp +++ b/flowexperimental/comp/flowexp_comp5.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<5>(int argc, char** argv) +int dispatchFlowExpComp<5, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp5_2p.cpp b/flowexperimental/comp/flowexp_comp5_2p.cpp new file mode 100644 index 000000000..c2c034d8a --- /dev/null +++ b/flowexperimental/comp/flowexp_comp5_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<5, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp6.cpp b/flowexperimental/comp/flowexp_comp6.cpp index b274591dd..95ee13066 100644 --- a/flowexperimental/comp/flowexp_comp6.cpp +++ b/flowexperimental/comp/flowexp_comp6.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<6>(int argc, char** argv) +int dispatchFlowExpComp<6, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp6_2p.cpp b/flowexperimental/comp/flowexp_comp6_2p.cpp new file mode 100644 index 000000000..3d6b88210 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp6_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<6, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp7.cpp b/flowexperimental/comp/flowexp_comp7.cpp index a40d2d23b..0c5fcc7f2 100644 --- a/flowexperimental/comp/flowexp_comp7.cpp +++ b/flowexperimental/comp/flowexp_comp7.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<7>(int argc, char** argv) +int dispatchFlowExpComp<7, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp7_2p.cpp b/flowexperimental/comp/flowexp_comp7_2p.cpp new file mode 100644 index 000000000..0dba55e3b --- /dev/null +++ b/flowexperimental/comp/flowexp_comp7_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, SINTEF Digital + + 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 "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<7, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/opm/models/blackoil/blackoilindices.hh b/opm/models/blackoil/blackoilindices.hh index 8350373d0..53f4ed6bb 100644 --- a/opm/models/blackoil/blackoilindices.hh +++ b/opm/models/blackoil/blackoilindices.hh @@ -94,10 +94,10 @@ struct BlackOilIndices numEnergy + numFoam + numBrine + numMICPs; //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx) + static constexpr int canonicalToActiveComponentIndex(const int compIdx) { return compIdx; } - static constexpr unsigned activeToCanonicalComponentIndex(unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex(const int compIdx) { return compIdx; } //////// diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index c2abedb02..863c8b00e 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -30,6 +30,8 @@ #include +#include + namespace Opm { /*! @@ -179,15 +181,15 @@ struct BlackOilOnePhaseIndices ////////////////////// //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(unsigned /*compIdx*/) + static constexpr int canonicalToActiveComponentIndex(const int /*compIdx*/) { return 0; } - static unsigned activeToCanonicalComponentIndex([[maybe_unused]] unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex([[maybe_unused]] const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; - assert(compIdx == 0); + constexpr_assert(compIdx == 0); if (gasEnabled) { return 2; } else if (waterEnabled) { diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index e4eefa50d..9d112a148 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/opm/models/blackoil/blackoiltwophaseindices.hh @@ -181,7 +181,7 @@ struct BlackOilTwoPhaseIndices ////////////////////// //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(const unsigned compIdx) + static constexpr int canonicalToActiveComponentIndex(const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; if (!gasEnabled) { @@ -201,10 +201,10 @@ struct BlackOilTwoPhaseIndices return compIdx - 1; } - static unsigned activeToCanonicalComponentIndex(unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex(const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; - assert(compIdx < 2); + constexpr_assert(compIdx < 2); if (!gasEnabled) { // oil = 0, water = 1 return compIdx; @@ -212,7 +212,7 @@ struct BlackOilTwoPhaseIndices // oil = 0, gas = 1 return compIdx * 2; } else { - assert(!oilEnabled); + constexpr_assert(!oilEnabled); } // water = 0, gas = 1; diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh index f067612f9..968a750ba 100644 --- a/opm/models/common/multiphasebaseproperties.hh +++ b/opm/models/common/multiphasebaseproperties.hh @@ -83,6 +83,8 @@ struct EnableDispersion { using type = UndefinedProperty; }; //! Enable convective mixing? template struct EnableConvectiveMixing { using type = UndefinedProperty; }; +template +struct EnableWater { using type = UndefinedProperty; }; } // namespace Opm::Properties diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh index 578b62d76..ac4c0ebf4 100644 --- a/opm/models/ptflash/flashindices.hh +++ b/opm/models/ptflash/flashindices.hh @@ -48,16 +48,20 @@ class FlashIndices { static constexpr int numComponents = getPropValue(); enum { enableEnergy = getPropValue() }; + enum { enableWater = getPropValue() }; using EnergyIndices = Opm::EnergyIndices; public: - static constexpr bool waterEnabled = false; + //! All phases active (note: immiscible/"dummy" water phase) + static constexpr bool waterEnabled = enableWater; static constexpr bool gasEnabled = true; static constexpr bool oilEnabled = true; - static constexpr int waterPhaseIdx = -1; - static constexpr int numPhases = 2; + + //! number of active phases + static constexpr int numPhases = enableWater ? 3 : 2; + //! number of equations/primary variables - static const int numEq = numComponents + EnergyIndices::numEq_; + static const int numEq = numComponents + EnergyIndices::numEq_ + (enableWater ? 1 : 0); // Primary variable indices @@ -66,6 +70,9 @@ public: //! Index of the molefraction of the first component static constexpr int z0Idx = pressure0Idx + 1; + + //! Index of water saturation + static constexpr int water0Idx = enableWater ? z0Idx + numComponents - 1 : -1000; // equation indices diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 7e71ef5a3..4237d19cb 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -76,6 +76,9 @@ class FlashIntensiveQuantities enum { enableEnergy = getPropValue() }; enum { dimWorld = GridView::dimensionworld }; enum { pressure0Idx = Indices::pressure0Idx }; + enum { water0Idx = Indices::water0Idx}; + + static constexpr bool waterEnabled = Indices::waterEnabled; using Scalar = GetPropType; using Evaluation = GetPropType; @@ -216,19 +219,26 @@ public: // Update saturation - // \Note: the current implementation assume oil-gas system. + Evaluation Sw = 0.0; + if constexpr (waterEnabled) { + Sw = priVars.makeEvaluation(water0Idx, timeIdx); + } Evaluation L = fluidState_.L(); - Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); - Evaluation Sg = Opm::max(1 - So, 0.0); - Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg); + Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); + Evaluation Sg = Opm::max(1 - So - Sw, 0.0); + Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw); So /= sumS; Sg /= sumS; - fluidState_.setSaturation(0, So); - fluidState_.setSaturation(1, Sg); + fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So); + fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg); + if constexpr (waterEnabled) { + Sw /= sumS; + fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); + } - fluidState_.setCompressFactor(0, Z_L); - fluidState_.setCompressFactor(1, Z_V); + fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L); + fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V); // Print saturation if (flashVerbosity >= 5) { @@ -250,7 +260,10 @@ public: // set the phase viscosity and density for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - paramCache.updatePhase(fluidState_, phaseIdx); + if (phaseIdx == static_cast(FluidSystem::oilPhaseIdx) + || phaseIdx == static_cast(FluidSystem::gasPhaseIdx)) { + paramCache.updatePhase(fluidState_, phaseIdx); + } const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx); diff --git a/opm/models/ptflash/flashlocalresidual.hh b/opm/models/ptflash/flashlocalresidual.hh index c528189e7..29ac021db 100644 --- a/opm/models/ptflash/flashlocalresidual.hh +++ b/opm/models/ptflash/flashlocalresidual.hh @@ -49,18 +49,24 @@ class FlashLocalResidual: public GetPropType; using IntensiveQuantities = GetPropType; using ElementContext = GetPropType; + using FluidSystem = GetPropType; enum { numEq = getPropValue() }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; + enum { water0Idx = Indices::water0Idx }; enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + enum { enableDiffusion = getPropValue() }; using DiffusionModule = Opm::DiffusionModule; enum { enableEnergy = getPropValue() }; using EnergyModule = Opm::EnergyModule; + static const bool waterEnabled = Indices::waterEnabled; + using Toolbox = Opm::MathToolbox; public: @@ -77,15 +83,25 @@ public: const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); const auto& fs = intQuants.fluidState(); - // compute storage term of all components within all phases - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - unsigned eqIdx = conti0EqIdx + compIdx; - storage[eqIdx] += - Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)) - * Toolbox::template decay(fs.density(phaseIdx)) + // compute water storage term + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { + unsigned eqIdx = conti0EqIdx + numComponents; + storage[eqIdx] = + Toolbox::template decay(fs.density(phaseIdx)) * Toolbox::template decay(fs.saturation(phaseIdx)) * Toolbox::template decay(intQuants.porosity()); } + else { + // compute storage term of all components within oil/gas phases + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + unsigned eqIdx = conti0EqIdx + compIdx; + storage[eqIdx] += + Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)) + * Toolbox::template decay(fs.density(phaseIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + } + } EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx); } @@ -146,19 +162,31 @@ public: up.fluidState().density(phaseIdx) * extQuants.volumeFlux(phaseIdx); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - flux[conti0EqIdx + compIdx] += - tmp*up.fluidState().massFraction(phaseIdx, compIdx); + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { + unsigned eqIdx = conti0EqIdx + numComponents; + flux[eqIdx] = tmp; + } + else { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*up.fluidState().massFraction(phaseIdx, compIdx); + } } } else { Evaluation tmp = Toolbox::value(up.fluidState().density(phaseIdx)) * extQuants.volumeFlux(phaseIdx); - - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - flux[conti0EqIdx + compIdx] += - tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx)); + + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { + unsigned eqIdx = conti0EqIdx + numComponents; + flux[eqIdx] = tmp; + } + else { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx)); + } } } } diff --git a/opm/models/ptflash/flashmodel.hh b/opm/models/ptflash/flashmodel.hh index f348523af..8a192f50e 100644 --- a/opm/models/ptflash/flashmodel.hh +++ b/opm/models/ptflash/flashmodel.hh @@ -135,6 +135,10 @@ template struct EnableEnergy { static constexpr bool value = false; }; +template +struct EnableWater +{ static constexpr int value = GetPropType::waterEnabled; }; + } // namespace Opm::Properties namespace Opm { diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh index 32500f38b..5b6fa1ddb 100644 --- a/opm/models/ptflash/flashnewtonmethod.hh +++ b/opm/models/ptflash/flashnewtonmethod.hh @@ -63,6 +63,8 @@ class FlashNewtonMethod : public GetPropType() }; + static constexpr bool waterEnabled = Indices::waterEnabled; + public: /*! * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& ) @@ -125,6 +127,14 @@ protected: for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol); } + + if constexpr (waterEnabled) { + // limit change in water saturation to 0.2 + constexpr Scalar dSwMax = 0.2; + if (update[Indices::water0Idx] > dSwMax) { + nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax; + } + } } private: void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh index 2df1c67c3..897c322ff 100644 --- a/opm/models/ptflash/flashprimaryvariables.hh +++ b/opm/models/ptflash/flashprimaryvariables.hh @@ -61,6 +61,14 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables // primary variable indices enum { z0Idx = Indices::z0Idx }; enum { pressure0Idx = Indices::pressure0Idx }; + enum { water0Idx = Indices::water0Idx }; + + static constexpr bool waterEnabled = Indices::waterEnabled; + + // phase indices + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; @@ -108,10 +116,17 @@ public: // the energy module EnergyModule::setPriVarTemperatures(*this, fluidState); + // assign components total fraction for (int i = 0; i < numComponents - 1; ++i) (*this)[z0Idx + i] = getValue(fluidState.moleFraction(i)); - (*this)[pressure0Idx] = getValue(fluidState.pressure(0)); + // assign pressure + (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx)); + + // assign water saturation + if constexpr (waterEnabled) { + (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx)); + } } /*! @@ -121,12 +136,15 @@ public: */ void print(std::ostream& os = std::cout) const { - os << "(p_" << FluidSystem::phaseName(0) << " = " + os << "(p_" << FluidSystem::phaseName(FluidSystem::oilPhaseIdx) << " = " << this->operator[](pressure0Idx); for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) { os << ", z_" << FluidSystem::componentName(compIdx) << " = " << this->operator[](z0Idx + compIdx); } + if constexpr (waterEnabled) { + os << ", S_w = " << this->operator[](water0Idx); + } os << ")" << std::flush; } }; diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp new file mode 100644 index 000000000..0a51ac903 --- /dev/null +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -0,0 +1,191 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include + +#include + +#include +#include + +#include + +namespace Opm { + +template +void CompositionalContainer:: +allocate(const unsigned bufferSize, + std::map& rstKeywords) +{ + if (auto& zmf = rstKeywords["ZMF"]; zmf > 0) { + this->allocated_ = true; + zmf = 0; + for (int i = 0; i < numComponents; ++i) { + moleFractions_[i].resize(bufferSize, 0.0); + } + } + + if (auto& xmf = rstKeywords["XMF"]; xmf > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) { + this->allocated_ = true; + xmf = 0; + for (int i = 0; i < numComponents; ++i) { + phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); + } + } + + if (auto& ymf = rstKeywords["YMF"]; ymf > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) { + this->allocated_ = true; + ymf = 0; + for (int i = 0; i < numComponents; ++i) { + phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); + } + } +} + +template +void CompositionalContainer:: +assignGasFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (phaseMoleFractions_[gasPhaseIdx][0].empty()) { + return; + } + + std::for_each(phaseMoleFractions_[gasPhaseIdx].begin(), + phaseMoleFractions_[gasPhaseIdx].end(), + [globalDofIdx, &fractions, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + +template +void CompositionalContainer:: +assignMoleFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (moleFractions_.empty()) { + return; + } + + std::for_each(moleFractions_.begin(), moleFractions_.end(), + [&fractions, globalDofIdx, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + +template +void CompositionalContainer:: +assignOilFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (phaseMoleFractions_[oilPhaseIdx][0].empty()) { + return; + } + + std::for_each(phaseMoleFractions_[oilPhaseIdx].begin(), + phaseMoleFractions_[oilPhaseIdx].end(), + [globalDofIdx, &fractions, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + +template +void CompositionalContainer:: +outputRestart(data::Solution& sol, + ScalarBuffer& oil_saturation) +{ + using DataEntry = + std::tuple&>; + + auto doInsert = [&sol](DataEntry& entry, + const data::TargetType target) + { + if (std::get<2>(entry).empty()) { + return; + } + + sol.insert(std::get(entry), + std::get(entry), + std::move(std::get<2>(entry)), + target); + }; + + auto entries = std::vector{}; + + // ZMF + if (!moleFractions_[0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]); + } + } + + // XMF + if (!phaseMoleFractions_[oilPhaseIdx][0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, + phaseMoleFractions_[oilPhaseIdx][i]); + } + } + + // YMF + if (!phaseMoleFractions_[gasPhaseIdx][0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, + phaseMoleFractions_[gasPhaseIdx][i]); + } + } + + if (!oil_saturation.empty()) { + entries.emplace_back("SOIL", UnitSystem::measure::identity, oil_saturation); + } + + std::for_each(entries.begin(), entries.end(), + [&doInsert](auto& array) + { doInsert(array, data::TargetType::RESTART_SOLUTION); }); + + this->allocated_ = false; +} + +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ + template class CompositionalContainer>; + +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class CompositionalContainer>; + +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP(2) +INSTANTIATE_COMP(3) +INSTANTIATE_COMP(4) +INSTANTIATE_COMP(5) +INSTANTIATE_COMP(6) +INSTANTIATE_COMP(7) + +} // namespace Opm diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp new file mode 100644 index 000000000..67a9ca4f9 --- /dev/null +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -0,0 +1,83 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::OutputBlackOilModule + */ +#ifndef OPM_COMPOSITIONAL_CONTAINER_HPP +#define OPM_COMPOSITIONAL_CONTAINER_HPP + +#include +#include +#include +#include +#include + +namespace Opm { + +namespace data { class Solution; } + +template +class CompositionalContainer +{ + using Scalar = typename FluidSystem::Scalar; + using ScalarBuffer = std::vector; + + static constexpr int numComponents = FluidSystem::numComponents; + + static constexpr int numPhases = FluidSystem::numPhases; + static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; + static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx; + static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; + +public: + void allocate(const unsigned bufferSize, + std::map& rstKeywords); + + using AssignFunction = std::function; + + void assignGasFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + + void assignMoleFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + + void assignOilFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + + void outputRestart(data::Solution& sol, + ScalarBuffer& oil_saturation); + + bool allocated() const + { return allocated_; } + +private: + bool allocated_ = false; + // total mole fractions for each component + std::array moleFractions_; + // mole fractions for each component in each phase + std::array, numPhases> phaseMoleFractions_; +}; + +} // namespace Opm + +#endif // OPM_COMPOSITIONAL_CONTAINER_HPP diff --git a/opm/simulators/flow/FIPContainer.cpp b/opm/simulators/flow/FIPContainer.cpp index 47bc1bb8f..95d6e37ef 100644 --- a/opm/simulators/flow/FIPContainer.cpp +++ b/opm/simulators/flow/FIPContainer.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include @@ -438,11 +438,19 @@ INSTANTIATE_TYPE(double) INSTANTIATE_TYPE(float) #endif -#define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasFluidSystem; \ +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class FIPContainer>; -INSTANTIATE_COMP(0) +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class FIPContainer>; + +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput INSTANTIATE_COMP(2) INSTANTIATE_COMP(3) INSTANTIATE_COMP(4) diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index be6a2616b..d6e6c817f 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -362,6 +362,9 @@ public: Dune::FieldVector z(0.0); Scalar sumMoles = 0.0; for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (Indices::waterEnabled && phaseIdx == static_cast(waterPhaseIdx)){ + continue; + } const auto saturation = fs.saturation(phaseIdx); for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation; @@ -491,7 +494,7 @@ protected: const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx); const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx); - if (water_active && Indices::numPhases > 1) + if (water_active && Indices::numPhases > 2) waterSaturationData = fp.get_double("SWAT"); else waterSaturationData.resize(numDof); @@ -527,6 +530,10 @@ protected: - waterSaturationData[dofIdx] - gasSaturationData[dofIdx]); } + if (water_active) { + dofFluidState.setSaturation(FluidSystem::waterPhaseIdx, + waterSaturationData[dofIdx]); + } ////// // set phase pressures diff --git a/opm/simulators/flow/FlowProblemCompProperties.hpp b/opm/simulators/flow/FlowProblemCompProperties.hpp index 12ecb2319..7b838a94c 100644 --- a/opm/simulators/flow/FlowProblemCompProperties.hpp +++ b/opm/simulators/flow/FlowProblemCompProperties.hpp @@ -69,17 +69,10 @@ private: using Scalar = GetPropType; using FluidSystem = GetPropType; - // using Traits = ThreePhaseMaterialTraits; - - // TODO: We should be able to use FluidSystem here and using Indices to handle the active phases - // some more development is needed using Traits = ThreePhaseMaterialTraits; + /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, + /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx, + /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>; public: using EclMaterialLawManager = ::Opm::EclMaterialLawManager; diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 46f556241..e7bef59a8 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -27,10 +27,10 @@ #include +#include #include #include - -#include +#include #include #include @@ -55,7 +55,6 @@ #include #include #include -#include #include #include #include @@ -104,8 +103,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState, bool enableBrine, bool enableSaltPrecipitation, bool enableExtbo, - bool enableMICP, - bool isCompositional) + bool enableMICP) : eclState_(eclState) , schedule_(schedule) , summaryState_(summaryState) @@ -124,7 +122,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState, , enableSaltPrecipitation_(enableSaltPrecipitation) , enableExtbo_(enableExtbo) , enableMICP_(enableMICP) - , isCompositional_(isCompositional) + , tracerC_(eclState_) , local_data_valid_(false) { const auto& fp = eclState_.fieldProps(); @@ -527,38 +525,6 @@ assignToSolution(data::Solution& sol) DataEntry{"TMULT_RC", UnitSystem::measure::identity, rockCompTransMultiplier_}, }; - // basically, for compositional, we can not use std::array for this. We need to generate the ZMF1, ZMF2, and so on - // and also, we need to map these values. - // TODO: the following should go to a function - if (this->isCompositional_) { - auto compositionalEntries = std::vector{}; - { - // ZMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]); - } - - // XMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - phaseMoleFractions_[oilPhaseIdx][i]); - } - - // YMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - phaseMoleFractions_[gasPhaseIdx][i]); - } - } - - for (auto& array: compositionalEntries) { - doInsert(array, data::TargetType::RESTART_SOLUTION); - } - } - for (auto& array : baseSolutionVector) { doInsert(array, data::TargetType::RESTART_SOLUTION); } @@ -602,15 +568,6 @@ assignToSolution(data::Solution& sol) data::TargetType::RESTART_SOLUTION); } - if (this->isCompositional_ && FluidSystem::phaseIsActive(oilPhaseIdx) && - ! this->saturation_[oilPhaseIdx].empty()) - { - sol.insert("SOIL", UnitSystem::measure::identity, - std::move(this->saturation_[oilPhaseIdx]), - data::TargetType::RESTART_SOLUTION); - } - - if ((eclState_.runspec().co2Storage() || eclState_.runspec().h2Storage()) && !rsw_.empty()) { auto mfrac = std::vector(this->rsw_.size(), 0.0); @@ -674,41 +631,7 @@ assignToSolution(data::Solution& sol) this->fipC_.outputRestart(sol); // Tracers - if (! this->freeTracerConcentrations_.empty()) { - const auto& tracers = this->eclState_.tracer(); - for (auto tracerIdx = 0*tracers.size(); - tracerIdx < tracers.size(); ++tracerIdx) - { - sol.insert(tracers[tracerIdx].fname(), - UnitSystem::measure::identity, - std::move(freeTracerConcentrations_[tracerIdx]), - data::TargetType::RESTART_TRACER_SOLUTION); - } - - // Put freeTracerConcentrations container into a valid state. Otherwise - // we'll move from vectors that have already been moved from if we - // get here and it's not a restart step. - this->freeTracerConcentrations_.clear(); - } - if (! this->solTracerConcentrations_.empty()) { - const auto& tracers = this->eclState_.tracer(); - for (auto tracerIdx = 0*tracers.size(); - tracerIdx < tracers.size(); ++tracerIdx) - { - if (solTracerConcentrations_[tracerIdx].empty()) - continue; - - sol.insert(tracers[tracerIdx].sname(), - UnitSystem::measure::identity, - std::move(solTracerConcentrations_[tracerIdx]), - data::TargetType::RESTART_TRACER_SOLUTION); - } - - // Put solTracerConcentrations container into a valid state. Otherwise - // we'll move from vectors that have already been moved from if we - // get here and it's not a restart step. - this->solTracerConcentrations_.clear(); - } + this->tracerC_.outputRestart(sol); } template @@ -828,16 +751,15 @@ doAllocBuffers(const unsigned bufferSize, const bool substep, const bool log, const bool isRestart, - const bool vapparsActive, - const bool enablePCHysteresis, - const bool enableNonWettingHysteresis, - const bool enableWettingHysteresis, - const unsigned numTracers, - const std::vector& enableSolTracers, - const unsigned numOutputNnc) + const EclHysteresisConfig* hysteresisConfig, + const unsigned numOutputNnc, + std::map rstKeywords) { + if (rstKeywords.empty()) { + rstKeywords = schedule_.rst_keywords(reportStepNum); + } + // Output RESTART_OPM_EXTENDED only when explicitly requested by user. - std::map rstKeywords = schedule_.rst_keywords(reportStepNum); for (auto& [keyword, should_write] : rstKeywords) { if (this->isOutputCreationDirective_(keyword)) { // 'BASIC', 'FREQ' and similar. Don't attempt to create @@ -1027,11 +949,13 @@ doAllocBuffers(const unsigned bufferSize, this->micpC_.allocate(bufferSize); } + const bool vapparsActive = schedule_[std::max(reportStepNum, 0u)].oilvap().getType() == + OilVaporizationProperties::OilVaporization::VAPPARS; if (vapparsActive) { soMax_.resize(bufferSize, 0.0); } - if (enableNonWettingHysteresis) { + if (hysteresisConfig && hysteresisConfig->enableNonWettingHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ soMax_.resize(bufferSize, 0.0); @@ -1043,7 +967,7 @@ doAllocBuffers(const unsigned bufferSize, //TODO add support for gas-water } } - if (enableWettingHysteresis) { + if (hysteresisConfig && hysteresisConfig->enableWettingHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ swMax_.resize(bufferSize, 0.0); @@ -1055,7 +979,7 @@ doAllocBuffers(const unsigned bufferSize, //TODO add support for gas-water } } - if (enablePCHysteresis) { + if (hysteresisConfig && hysteresisConfig->enablePCHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ swmin_.resize(bufferSize, 0.0); @@ -1259,19 +1183,7 @@ doAllocBuffers(const unsigned bufferSize, } // tracers - if (numTracers > 0) { - freeTracerConcentrations_.resize(numTracers); - for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) - { - freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); - } - solTracerConcentrations_.resize(numTracers); - for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) - { - if (enableSolTracers[tracerIdx]) - solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); - } - } + this->tracerC_.allocate(bufferSize); if (rstKeywords["RESIDUAL"] > 0) { rstKeywords["RESIDUAL"] = 0; @@ -1293,30 +1205,6 @@ doAllocBuffers(const unsigned bufferSize, overburdenPressure_.resize(bufferSize, 0.0); } - if (this->isCompositional_) { - if (rstKeywords["ZMF"] > 0) { - rstKeywords["ZMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - moleFractions_[i].resize(bufferSize, 0.0); - } - } - - if (rstKeywords["XMF"] > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) { - rstKeywords["XMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); - } - } - - if (rstKeywords["YMF"] > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) { - rstKeywords["YMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); - } - } - } - - //Warn for any unhandled keyword if (log) { for (auto& keyValue: rstKeywords) { @@ -1585,12 +1473,19 @@ INSTANTIATE_TYPE(double) INSTANTIATE_TYPE(float) #endif -#define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasFluidSystem; \ +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class GenericOutputBlackoilModule>; -INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class GenericOutputBlackoilModule>; +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput INSTANTIATE_COMP(2) INSTANTIATE_COMP(3) INSTANTIATE_COMP(4) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index c2d45aafa..b73cc4144 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include @@ -62,6 +63,7 @@ struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false; namespace Opm { namespace data { class Solution; } +class EclHysteresisConfig; class EclipseState; class Schedule; class SummaryConfig; @@ -319,21 +321,16 @@ protected: bool enableBrine, bool enableSaltPrecipitation, bool enableExtbo, - bool enableMICP, - bool isCompositional = false); + bool enableMICP); void doAllocBuffers(unsigned bufferSize, unsigned reportStepNum, const bool substep, const bool log, const bool isRestart, - const bool vapparsActive = false, - const bool enablePCHysteresis = false, - const bool enableNonWettingHysteresis =false, - const bool enableWettingHysteresis = false, - unsigned numTracers = 0, - const std::vector& enableSolTracers = {}, - unsigned numOutputNnc = 0); + const EclHysteresisConfig* hysteresisConfig, + unsigned numOutputNnc = 0, + std::map rstKeywords = {}); void makeRegionSum(Inplace& inplace, const std::string& region_name, @@ -390,7 +387,6 @@ protected: bool enableSaltPrecipitation_{false}; bool enableExtbo_{false}; bool enableMICP_{false}; - bool isCompositional_{false}; bool forceDisableFipOutput_{false}; bool forceDisableFipresvOutput_{false}; @@ -468,12 +464,7 @@ protected: std::array viscosity_; std::array relativePermeability_; - // total mole fractions for each component - std::array moleFractions_; - // mole fractions for each component in each phase - std::array, numPhases> phaseMoleFractions_; - std::vector freeTracerConcentrations_; - std::vector solTracerConcentrations_; + TracerContainer tracerC_; std::array residual_; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 650cb8f3e..a6fcc58f2 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -184,12 +184,7 @@ public: substep, log, isRestart, - problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)), - problem.materialLawManager()->enablePCHysteresis(), - problem.materialLawManager()->enableNonWettingHysteresis(), - problem.materialLawManager()->enableWettingHysteresis(), - problem.tracerModel().numTracers(), - problem.tracerModel().enableSolTracers(), + &problem.materialLawManager()->hysteresisConfig(), problem.eclWriter()->getOutputNnc().size()); } @@ -517,7 +512,7 @@ public: const auto& matLawManager = problem.materialLawManager(); if (matLawManager->enableHysteresis()) { - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) { Scalar somax; Scalar swmax; @@ -525,7 +520,7 @@ public: matLawManager->oilWaterHysteresisParams( somax, swmax, swmin, globalDofIdx); - + if (matLawManager->enableNonWettingHysteresis()) { if (!this->soMax_.empty()) { this->soMax_[globalDofIdx] = somax; @@ -543,14 +538,14 @@ public: } } - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { Scalar sgmax; Scalar shmax; Scalar somin; matLawManager->gasOilHysteresisParams( sgmax, shmax, somin, globalDofIdx); - + if (matLawManager->enableNonWettingHysteresis()) { if (!this->sgmax_.empty()) { this->sgmax_[globalDofIdx] = sgmax; @@ -568,7 +563,7 @@ public: } } } else { - + if (!this->soMax_.empty()) this->soMax_[globalDofIdx] = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx)); @@ -647,25 +642,14 @@ public: // tracers const auto& tracerModel = simulator_.problem().tracerModel(); - if (! this->freeTracerConcentrations_.empty()) { - for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->freeTracerConcentrations_[tracerIdx].empty()) { - continue; - } - this->freeTracerConcentrations_[tracerIdx][globalDofIdx] = - tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); - } - } - if (! this->solTracerConcentrations_.empty()) { - for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->solTracerConcentrations_[tracerIdx].empty()) { - continue; - } - this->solTracerConcentrations_[tracerIdx][globalDofIdx] = - tracerModel.solTracerConcentration(tracerIdx, globalDofIdx); - - } - } + this->tracerC_.assignFreeConcentrations(globalDofIdx, + [globalDofIdx, &tracerModel](const unsigned tracerIdx) + { return tracerModel.freeTracerConcentration(tracerIdx, + globalDofIdx); }); + this->tracerC_.assignSolConcentrations(globalDofIdx, + [globalDofIdx, &tracerModel](const unsigned tracerIdx) + { return tracerModel.solTracerConcentration(tracerIdx, + globalDofIdx); }); // output residual for ( int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx ) @@ -1176,7 +1160,7 @@ public: if (simulator.problem().materialLawManager()->enableHysteresis()) { auto matLawManager = simulator.problem().materialLawManager(); - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) { Scalar somax = 2.0; Scalar swmax = -2.0; @@ -1200,7 +1184,7 @@ public: matLawManager->setOilWaterHysteresisParams( somax, swmax, swmin, elemIdx); } - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { Scalar sgmax = 2.0; Scalar shmax = -2.0; @@ -1538,7 +1522,7 @@ private: { this->updateCO2InGas(globalDofIdx, pv, intQuants); } - + if (this->fipC_.hasCo2InWater() && (FluidSystem::phaseIsActive(waterPhaseIdx) || FluidSystem::phaseIsActive(oilPhaseIdx))) diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index d6102528d..97606d9d8 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -32,15 +32,20 @@ #include #include +#include #include #include #include #include + +#include +#include #include #include +#include #include #include @@ -53,8 +58,7 @@ #include -namespace Opm -{ +namespace Opm { // forward declaration template @@ -102,8 +106,7 @@ public: getPropValue(), getPropValue(), getPropValue(), - getPropValue(), - true) + getPropValue()) , simulator_(simulator) { for (auto& region_pair : this->regions_) { @@ -155,11 +158,19 @@ public: return; } - this->doAllocBuffers(bufferSize, - reportStepNum, - substep, - log, - isRestart); + auto rstKeywords = this->schedule_.rst_keywords(reportStepNum); + this->compC_.allocate(bufferSize, rstKeywords); + + this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart, + /* hysteresisConfig = */ nullptr, + /* numOutputNnc =*/ 0, + std::move(rstKeywords)); + } + + void assignToSolution(data::Solution& sol) + { + this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]); + BaseType::assignToSolution(sol); } /*! @@ -187,20 +198,21 @@ public: Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]); } - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (this->moleFractions_[compIdx].empty()) continue; + if (this->compC_.allocated()) { + this->compC_.assignMoleFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(compIdx)); }); - this->moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx)); - } - // XMF and YMF - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - if (this->phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue; - this->phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx)); - } if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - if (this->phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue; - this->phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx)); + this->compC_.assignGasFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(gasPhaseIdx, compIdx)); }); + } + + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + this->compC_.assignOilFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(oilPhaseIdx, compIdx)); }); } } @@ -332,6 +344,7 @@ private: } const Simulator& simulator_; + CompositionalContainer compC_; }; } // namespace Opm diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp new file mode 100644 index 000000000..9f73e758c --- /dev/null +++ b/opm/simulators/flow/TracerContainer.cpp @@ -0,0 +1,154 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include + +namespace Opm { + +template +void TracerContainer:: +allocate(const unsigned bufferSize) +{ + const auto& tracers = eclState_.tracer(); + if (!tracers.empty()) { + allocated_ = true; + freeConcentrations_.resize(tracers.size()); + solConcentrations_.resize(tracers.size()); + std::for_each(tracers.begin(), tracers.end(), + [idx = 0, bufferSize, this](const auto& tracer) mutable + { + freeConcentrations_[idx].resize(bufferSize, 0.0); + if (((tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) || + (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil())) && + (tracer.solution_concentration.has_value() || + tracer.solution_tvdp.has_value())) + { + solConcentrations_[idx].resize(bufferSize, 0.0); + } + ++idx; + }); + } +} + +template +void TracerContainer:: +assignFreeConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration) +{ + std::for_each(freeConcentrations_.begin(), freeConcentrations_.end(), + [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable + { + if (!tracer.empty()) { + tracer[globalDofIdx] = concentration(idx); + } + ++idx; + }); +} + +template +void TracerContainer:: +assignSolConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration) +{ + std::for_each(solConcentrations_.begin(), solConcentrations_.end(), + [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable + { + if (!tracer.empty()) { + tracer[globalDofIdx] = concentration(idx); + } + ++idx; + }); +} + +template +void TracerContainer:: +outputRestart(data::Solution& sol) +{ + if (!this->allocated_) { + return; + } + + const auto& tracers = this->eclState_.tracer(); + std::for_each(tracers.begin(), tracers.end(), + [idx = 0, &sol, this](const auto& tracer) mutable + { + sol.insert(tracer.fname(), + UnitSystem::measure::identity, + std::move(freeConcentrations_[idx]), + data::TargetType::RESTART_TRACER_SOLUTION); + + if (!solConcentrations_[idx].empty()) { + sol.insert(tracer.sname(), + UnitSystem::measure::identity, + std::move(solConcentrations_[idx]), + data::TargetType::RESTART_TRACER_SOLUTION); + } + ++idx; + }); + + this->allocated_ = false; +} + +template using FS = BlackOilFluidSystem; + +#define INSTANTIATE_TYPE(T) \ + template class TracerContainer>; + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ + template class TracerContainer>; + +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class TracerContainer>; + +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP(2) +INSTANTIATE_COMP(3) +INSTANTIATE_COMP(4) +INSTANTIATE_COMP(5) +INSTANTIATE_COMP(6) +INSTANTIATE_COMP(7) + +} // namespace Opm diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp new file mode 100644 index 000000000..dfd8ff86a --- /dev/null +++ b/opm/simulators/flow/TracerContainer.hpp @@ -0,0 +1,70 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::OutputBlackOilModule + */ +#ifndef OPM_TRACER_CONTAINER_HPP +#define OPM_TRACER_CONTAINER_HPP + +#include +#include + +namespace Opm { + +namespace data { class Solution; } +class EclipseState; + +template +class TracerContainer +{ + using Scalar = typename FluidSystem::Scalar; + using ScalarBuffer = std::vector; + +public: + TracerContainer(const EclipseState& eclState) + : eclState_(eclState) + {} + + void allocate(const unsigned bufferSize); + + using AssignFunction = std::function; + + void assignFreeConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration); + + void assignSolConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration); + + void outputRestart(data::Solution& sol); + +private: + const EclipseState& eclState_; + + std::vector freeConcentrations_{}; + std::vector solConcentrations_{}; + bool allocated_{false}; +}; + +} // namespace Opm + +#endif // OPM_TRACER_CONTAINER_HPP diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 053a6fbc6..e118fb647 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -377,10 +378,11 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverPrepare); + try { + initPrepare(M,b); - initPrepare(M,b); - - prepareFlexibleSolver(); + prepareFlexibleSolver(); + } OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. Check for errors related to missing nodes."); } diff --git a/opm/simulators/linalg/gpuistl/GpuBuffer.cpp b/opm/simulators/linalg/gpuistl/GpuBuffer.cpp index b4132ce09..49ce121e1 100644 --- a/opm/simulators/linalg/gpuistl/GpuBuffer.cpp +++ b/opm/simulators/linalg/gpuistl/GpuBuffer.cpp @@ -188,6 +188,7 @@ GpuBuffer::copyFromHost(const std::vector& data) { copyFromHost(data.data(), data.size()); } + template void GpuBuffer::copyToHost(std::vector& data) const diff --git a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp index 369ea3ba2..1ba5cb861 100644 --- a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp +++ b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp @@ -27,6 +27,7 @@ #include #include #include +#include namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.cpp b/opm/simulators/linalg/gpuistl/GpuDILU.cpp index d311386b0..2011fb2a2 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.cpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.cpp @@ -37,6 +37,7 @@ #include #include #include + namespace Opm::gpuistl { @@ -93,7 +94,8 @@ GpuDILU::GpuDILU(const M& A, bool splitMatrix, bool tuneKernels, int } } - computeDiagAndMoveReorderedData(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize); + reorderAndSplitMatrix(m_moveThreadBlockSize); + computeDiagonal(m_DILUFactorizationThreadBlockSize); if (m_tuneThreadBlockSizes) { tuneThreadBlockSizes(); @@ -110,10 +112,31 @@ template void GpuDILU::apply(X& v, const Y& d) { + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + OPM_TIMEBLOCK(prec_apply); { - apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + const auto ptrs = std::make_pair(v.data(), d.data()); + + auto it = m_apply_graphs.find(ptrs); + + if (it == m_apply_graphs.end()) { + OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal)); + + // The apply functions contains lots of small function calls which call a kernel each + apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + + OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get())); + OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); + } + OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0)); } + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } template @@ -135,7 +158,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInvFloat->data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { detail::DILU::solveLowerLevelSetSplit( m_gpuMatrixReorderedLowerFloat->getNonZeroValues().data(), @@ -147,7 +171,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else { detail::DILU::solveLowerLevelSetSplit( m_gpuMatrixReorderedLower->getNonZeroValues().data(), @@ -159,7 +184,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } } else { detail::DILU::solveLowerLevelSet( @@ -172,7 +198,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } levelStartIdx += numOfRowsInLevel; } @@ -193,7 +220,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInvFloat->data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG){ detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(), @@ -204,7 +232,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else { detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpper->getNonZeroValues().data(), @@ -215,7 +244,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } else { detail::DILU::solveUpperLevelSet( @@ -227,7 +257,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } } @@ -252,7 +283,9 @@ GpuDILU::update() { OPM_TIMEBLOCK(prec_update); { - update(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize); + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(m_moveThreadBlockSize); + computeDiagonal(m_DILUFactorizationThreadBlockSize); } } @@ -260,13 +293,22 @@ template void GpuDILU::update(int moveThreadBlockSize, int factorizationBlockSize) { - m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu - computeDiagAndMoveReorderedData(moveThreadBlockSize, factorizationBlockSize); + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(moveThreadBlockSize); + computeDiagonal(factorizationBlockSize); + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } template void -GpuDILU::computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationBlockSize) +GpuDILU::reorderAndSplitMatrix(int moveThreadBlockSize) { if (m_splitMatrix) { detail::copyMatDataToReorderedSplit( @@ -290,7 +332,12 @@ GpuDILU::computeDiagAndMoveReorderedData(int moveThreadBlockSize, in m_gpuMatrixReordered->N(), moveThreadBlockSize); } +} +template +void +GpuDILU::computeDiagonal(int factorizationBlockSize) +{ int levelStartIdx = 0; for (int level = 0; level < m_levelSets.size(); ++level) { const int numOfRowsInLevel = m_levelSets[level].size(); diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.hpp b/opm/simulators/linalg/gpuistl/GpuDILU.hpp index 241eef3e9..a023c05d4 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.hpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.hpp @@ -24,8 +24,10 @@ #include #include #include +#include #include - +#include +#include namespace Opm::gpuistl @@ -82,8 +84,11 @@ public: //! \brief Updates the matrix data. void update() final; + //! \brief perform matrix splitting and reordering + void reorderAndSplitMatrix(int moveThreadBlockSize); + //! \brief Compute the diagonal of the DILU, and update the data of the reordered matrix - void computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationThreadBlockSize); + void computeDiagonal(int factorizationThreadBlockSize); //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels void tuneThreadBlockSizes(); @@ -153,6 +158,16 @@ private: int m_lowerSolveThreadBlockSize = -1; int m_moveThreadBlockSize = -1; int m_DILUFactorizationThreadBlockSize = -1; + + // Graphs for Apply + std::map, GPUGraph> m_apply_graphs; + std::map, GPUGraphExec> m_executableGraphs; + + // Stream for the DILU operations on the GPU + GPUStream m_stream{}; + // Events for synchronization with main stream + GPUEvent m_before{}; + GPUEvent m_after{}; }; } // end namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/GpuVector.cpp b/opm/simulators/linalg/gpuistl/GpuVector.cpp index ddafc1dbf..af39e9941 100644 --- a/opm/simulators/linalg/gpuistl/GpuVector.cpp +++ b/opm/simulators/linalg/gpuistl/GpuVector.cpp @@ -266,6 +266,19 @@ GpuVector::copyFromHost(const T* dataPointer, size_t numberOfElements) OPM_GPU_SAFE_CALL(cudaMemcpy(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice)); } +template +void +GpuVector::copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream) +{ + if (numberOfElements > dim()) { + OPM_THROW(std::runtime_error, + fmt::format("Requesting to copy too many elements. Vector has {} elements, while {} was requested.", + dim(), + numberOfElements)); + } + OPM_GPU_SAFE_CALL(cudaMemcpyAsync(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice, stream)); +} + template void GpuVector::copyToHost(T* dataPointer, size_t numberOfElements) const diff --git a/opm/simulators/linalg/gpuistl/GpuVector.hpp b/opm/simulators/linalg/gpuistl/GpuVector.hpp index e369d1d9d..66a1ebd65 100644 --- a/opm/simulators/linalg/gpuistl/GpuVector.hpp +++ b/opm/simulators/linalg/gpuistl/GpuVector.hpp @@ -203,6 +203,7 @@ public: * @note assumes that this vector has numberOfElements elements */ void copyFromHost(const T* dataPointer, size_t numberOfElements); + void copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream); /** * @brief copyFromHost copies numberOfElements to the CPU memory dataPointer diff --git a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp index 686b8f09f..42b535b76 100644 --- a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp +++ b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp @@ -37,6 +37,7 @@ #include #include #include + namespace Opm::gpuistl { @@ -58,7 +59,6 @@ OpmGpuILU0::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel , m_tuneThreadBlockSizes(tuneKernels) , m_mixedPrecisionScheme(makeMatrixStorageMPScheme(mixedPrecisionScheme)) { - // TODO: Should in some way verify that this matrix is symmetric, only do it debug mode? // Some sanity check OPM_ERROR_IF(A.N() != m_gpuMatrix.N(), @@ -98,7 +98,8 @@ OpmGpuILU0::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel } } - LUFactorizeAndMoveData(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize); + reorderAndSplitMatrix(m_moveThreadBlockSize); + LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize); if (m_tuneThreadBlockSizes) { tuneThreadBlockSizes(); @@ -117,7 +118,29 @@ OpmGpuILU0::apply(X& v, const Y& d) { OPM_TIMEBLOCK(prec_apply); { - apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + + const auto ptrs = std::make_pair(v.data(), d.data()); + + auto it = m_apply_graphs.find(ptrs); + + if (it == m_apply_graphs.end()) { + OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal)); + + // The apply functions contains lots of small function calls which call a kernel each + apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + + OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get())); + OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); + } + OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0)); + + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } } @@ -142,7 +165,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else{ detail::ILU0::solveLowerLevelSetSplit( @@ -154,7 +178,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } } else { detail::ILU0::solveLowerLevelSet(m_gpuReorderedLU->getNonZeroValues().data(), @@ -165,7 +190,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } levelStartIdx += numOfRowsInLevel; } @@ -185,7 +211,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiagFloat.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { detail::ILU0::solveUpperLevelSetSplit( @@ -197,7 +224,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiag.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else{ detail::ILU0::solveUpperLevelSetSplit( @@ -209,7 +237,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiag.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } else { detail::ILU0::solveUpperLevelSet(m_gpuReorderedLU->getNonZeroValues().data(), @@ -219,7 +248,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i levelStartIdx, numOfRowsInLevel, v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } } @@ -241,10 +271,9 @@ template void OpmGpuILU0::update() { - OPM_TIMEBLOCK(prec_update); - { - update(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize); - } + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(m_moveThreadBlockSize); + LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize); } template @@ -253,13 +282,23 @@ OpmGpuILU0::update(int moveThreadBlockSize, int factorizationThreadB { OPM_TIMEBLOCK(prec_update); { + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu - LUFactorizeAndMoveData(moveThreadBlockSize, factorizationThreadBlockSize); + reorderAndSplitMatrix(moveThreadBlockSize); + LUFactorizeMatrix(factorizationThreadBlockSize); + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } } + template void -OpmGpuILU0::LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize) +OpmGpuILU0::reorderAndSplitMatrix(int moveThreadBlockSize) { if (m_splitMatrix) { detail::copyMatDataToReorderedSplit( @@ -283,6 +322,12 @@ OpmGpuILU0::LUFactorizeAndMoveData(int moveThreadBlockSize, int fact m_gpuReorderedLU->N(), moveThreadBlockSize); } +} + +template +void +OpmGpuILU0::LUFactorizeMatrix(int factorizationThreadBlockSize) +{ int levelStartIdx = 0; for (int level = 0; level < m_levelSets.size(); ++level) { const int numOfRowsInLevel = m_levelSets[level].size(); diff --git a/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp b/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp index c163e0483..e4cc3f059 100644 --- a/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp +++ b/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -84,8 +85,11 @@ public: //! \brief Updates the matrix data. void update() final; + //! \brief perform matrix splitting and reordering + void reorderAndSplitMatrix(int moveThreadBlockSize); + //! \brief Compute LU factorization, and update the data of the reordered matrix - void LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize); + void LUFactorizeMatrix(int factorizationThreadBlockSize); //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels void tuneThreadBlockSizes(); @@ -152,6 +156,16 @@ private: int m_lowerSolveThreadBlockSize = -1; int m_moveThreadBlockSize = -1; int m_ILU0FactorizationThreadBlockSize = -1; + + // Graphs for Apply + std::map, GPUGraph> m_apply_graphs; + std::map, GPUGraphExec> m_executableGraphs; + + // Stream for the DILU operations on the GPU + GPUStream m_stream{}; + // Events for synchronization with main stream + GPUEvent m_before{}; + GPUEvent m_after{}; }; } // end namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu index b8f9f48ca..8af6765df 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu @@ -85,10 +85,12 @@ namespace // TODO: removce the first condition in the for loop for (int block = nnzIdx; block < nnzIdxLim; ++block) { const int col = colIndices[block]; - mmvMixedGeneral(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs); + mmvMixedGeneral( + &mat[block * blocksize * blocksize], &v[col * blocksize], rhs); } - mvMixedGeneral(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); + mvMixedGeneral( + &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); } } @@ -137,10 +139,12 @@ namespace LinearSolverScalar rhs[blocksize] = {0}; for (int block = nnzIdx; block < nnzIdxLim; ++block) { const int col = colIndices[block]; - umvMixedGeneral(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs); + umvMixedGeneral( + &mat[block * blocksize * blocksize], &v[col * blocksize], rhs); } - mmvMixedGeneral(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); + mmvMixedGeneral( + &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); } } @@ -211,8 +215,8 @@ namespace } } - // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate results - // TOOD: The important part is to only cast after that is fully computed + // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate + // results TOOD: The important part is to only cast after that is fully computed template __global__ void cuComputeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, int* lowerRowIndices, @@ -239,7 +243,8 @@ namespace InputScalar dInvTmp[blocksize * blocksize]; for (int i = 0; i < blocksize; ++i) { for (int j = 0; j < blocksize; ++j) { - dInvTmp[i * blocksize + j] = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j]; + dInvTmp[i * blocksize + j] + = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j]; } } @@ -257,21 +262,26 @@ namespace if constexpr (detail::storeOffDiagonalAsFloat(mixedPrecisionScheme)) { // TODO: think long and hard about whether this performs only the wanted memory transfers - moveBlock(&srcReorderedLowerMat[block * blocksize * blocksize], &dstLowerMat[block * blocksize * blocksize]); - moveBlock(&srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], &dstUpperMat[symOppositeBlock * blocksize * blocksize]); + moveBlock( + &srcReorderedLowerMat[block * blocksize * blocksize], + &dstLowerMat[block * blocksize * blocksize]); + moveBlock( + &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], + &dstUpperMat[symOppositeBlock * blocksize * blocksize]); } mmx2Subtraction(&srcReorderedLowerMat[block * blocksize * blocksize], - &dInv[col * blocksize * blocksize], - &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], - dInvTmp); + &dInv[col * blocksize * blocksize], + &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], + dInvTmp); } invBlockInPlace(dInvTmp); moveBlock(dInvTmp, &dInv[reorderedRowIdx * blocksize * blocksize]); if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) { - moveBlock(dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important! + moveBlock( + dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important! } } } @@ -289,12 +299,13 @@ solveLowerLevelSet(T* reorderedMat, const T* dInv, const T* d, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSet<<>>( + cuSolveLowerLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); } @@ -310,13 +321,15 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat, const DiagonalScalar* dInv, const LinearSolverScalar* d, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveLowerLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); + cuSolveLowerLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); } // perform the upper solve for all rows in the same level set template @@ -329,12 +342,13 @@ solveUpperLevelSet(T* reorderedMat, int rowsInLevelSet, const T* dInv, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSet<<>>( + cuSolveUpperLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } @@ -348,13 +362,15 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveUpperLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); + cuSolveUpperLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } template @@ -409,43 +425,134 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuComputeDiluDiagonalSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuComputeDiluDiagonalSplit<<>>(srcReorderedLowerMat, - lowerRowIndices, - lowerColIndices, - srcReorderedUpperMat, - upperRowIndices, - upperColIndices, - srcDiagonal, - reorderedToNatural, - naturalToReordered, - startIdx, - rowsInLevelSet, - dInv, - dstDiag, - dstLowerMat, - dstUpperMat); + cuComputeDiluDiagonalSplit + <<>>(srcReorderedLowerMat, + lowerRowIndices, + lowerColIndices, + srcReorderedUpperMat, + upperRowIndices, + upperColIndices, + srcDiagonal, + reorderedToNatural, + naturalToReordered, + startIdx, + rowsInLevelSet, + dInv, + dstDiag, + dstLowerMat, + dstUpperMat); } else { OPM_THROW(std::invalid_argument, "Inverting diagonal is not implemented for blocksizes > 3"); } } -// TODO: format #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize) \ template void computeDiluDiagonal(T*, int*, int*, int*, int*, const int, int, T*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void solveUpperLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int); \ - template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, const T*, T*, int); + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void solveUpperLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \ + template void solveLowerLevelSet( \ + T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); INSTANTIATE_KERNEL_WRAPPERS(float, 1); INSTANTIATE_KERNEL_WRAPPERS(float, 2); @@ -460,19 +567,29 @@ INSTANTIATE_KERNEL_WRAPPERS(double, 4); INSTANTIATE_KERNEL_WRAPPERS(double, 5); INSTANTIATE_KERNEL_WRAPPERS(double, 6); -#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar) \ - template void solveUpperLevelSetSplit( \ - MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int); \ - template void solveLowerLevelSetSplit( \ - MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, const LinearSolverScalar*, LinearSolverScalar*, int); +#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar) \ + template void solveUpperLevelSetSplit( \ + MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int, cudaStream_t); \ + template void solveLowerLevelSetSplit( \ + MatrixScalar*, \ + int*, \ + int*, \ + int*, \ + int, \ + int, \ + const DiagonalScalar*, \ + const LinearSolverScalar*, \ + LinearSolverScalar*, \ + int, \ + cudaStream_t); // TODO: be smarter about this... Surely this instantiates many more combinations that are actually needed -#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(blocksize) \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, double); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, double); \ +#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(blocksize) \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, double); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, double); \ INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, double); INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(1); diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp index b7141834f..77a9cae69 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp @@ -54,7 +54,8 @@ void solveLowerLevelSet(T* reorderedMat, const T* dInv, const T* d, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel @@ -82,7 +83,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedUpperMat, const DiagonalScalar* dInv, const LinearSolverScalar* d, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -108,7 +110,8 @@ void solveUpperLevelSet(T* reorderedMat, int rowsInLevelSet, const T* dInv, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -134,7 +137,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedUpperMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Computes the ILU0 of the diagonal elements of the reordered matrix and stores it in a reordered vector diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu index 46183efd2..e66c42f8a 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu @@ -215,9 +215,9 @@ namespace // as of now, if we are using mixed precision, then we are always storing the off-diagonals as floats, // and sometimes also the diagonal. if constexpr (detail::usingMixedPrecision(mixedPrecisionScheme)) { - // if we are want to store the entire matrix as a float then we must also move the diagonal block from double to float - // if not then we just use the double diagonal that is already now stored in srcDiagonal - if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)){ + // if we are want to store the entire matrix as a float then we must also move the diagonal block from + // double to float if not then we just use the double diagonal that is already now stored in srcDiagonal + if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) { moveBlock(&srcDiagonal[reorderedIdx * scalarsInBlock], &dstDiagonal[reorderedIdx * scalarsInBlock]); } @@ -362,12 +362,13 @@ solveLowerLevelSet(T* reorderedMat, int rowsInLevelSet, const T* d, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSet<<>>( + cuSolveLowerLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); } // perform the upper solve for all rows in the same level set @@ -380,12 +381,13 @@ solveUpperLevelSet(T* reorderedMat, int startIdx, int rowsInLevelSet, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSet<<>>( + cuSolveUpperLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, v); } @@ -399,13 +401,15 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const LinearSolverScalar* d, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveLowerLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); + cuSolveLowerLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); } // perform the upper solve for all rows in the same level set template @@ -418,13 +422,15 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveUpperLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); + cuSolveUpperLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } template @@ -484,28 +490,28 @@ LUFactorizationSplit(InputScalar* srcReorderedLowerMat, } #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize) \ - template void solveUpperLevelSet(T*, int*, int*, int*, int, int, T*, int); \ - template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int); \ + template void solveUpperLevelSet(T*, int*, int*, int*, int, int, T*, int, cudaStream_t); \ + template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \ template void LUFactorization(T*, int*, int*, int*, int*, size_t, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ - T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ + T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ - T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); + template void LUFactorizationSplit( \ + T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); -#define INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(T) \ - INSTANTIATE_KERNEL_WRAPPERS(T, 1); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 2); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 3); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 4); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 5); \ +#define INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(T) \ + INSTANTIATE_KERNEL_WRAPPERS(T, 1); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 2); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 3); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 4); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 5); \ INSTANTIATE_KERNEL_WRAPPERS(T, 6); INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(float) @@ -514,32 +520,32 @@ INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(double) #define INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(blocksize) \ /* double preconditioner */ \ template void solveLowerLevelSetSplit( \ - double*, int*, int*, int*, int, int, const double*, double*, int); \ + double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ template void solveLowerLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, double*, int); \ + float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float preconditioner */ \ template void solveLowerLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, float*, int); \ + float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t); \ \ /* double preconditioner */ \ - template void solveUpperLevelSetSplit( \ - double*, int*, int*, int*, int, int, const double*, double*, int); \ + template void solveUpperLevelSetSplit( \ + double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, double*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, float*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const double*, float*, int, cudaStream_t); \ /* double preconditioner */ \ - template void solveUpperLevelSetSplit( \ - double*, int*, int*, int*, int, int, const float*, double*, int); \ + template void solveUpperLevelSetSplit( \ + double*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, double*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t); \ /* float preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, float*, int); + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t); INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(1); diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp index cdd78ec91..04390ff77 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp @@ -46,7 +46,8 @@ void solveUpperLevelSet(T* reorderedMat, int startIdx, int rowsInLevelSet, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel @@ -72,7 +73,8 @@ void solveLowerLevelSet(T* reorderedMat, int rowsInLevelSet, const T* d, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -99,7 +101,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an lower solve on certain rows in a matrix that can safely be computed in parallel @@ -127,7 +130,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedLowerMat, int rowsInLevelSet, const LinearSolverScalar* d, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Computes the ILU Factorization of the input bcsr matrix, which is stored in a reordered way. The diagonal diff --git a/opm/simulators/utils/FullySupportedFlowKeywords.cpp b/opm/simulators/utils/FullySupportedFlowKeywords.cpp index ea13034b6..669635894 100644 --- a/opm/simulators/utils/FullySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/FullySupportedFlowKeywords.cpp @@ -32,6 +32,12 @@ const SupportedKeywordItems& fullySupported() { static const SupportedKeywordItems fully_supported_keywords_strings = { + { + "GEFAC", + { + {3,{true, is_bool_convertible {}, "GEFAC(GRPNETWK): String value must be convertible to bool."}}, // USE_GEFAC_IN_NETWORK + }, + }, { "NEXTSTEP", { @@ -44,6 +50,12 @@ fullySupported() {3,{true, allow_values {"ORAT", "WRAT", "GRAT", "LRAT", "RESV", "BHP"}, "WCONHIST(TARGET): should be set to ORAT/WRAT/GRAT/LRAT/RESV or BHP"}}, // CMODE }, }, + { + "WEFAC", + { + {3,{true, is_bool_convertible {}, "WEFAC(WELNETWK): String value must be convertible to bool."}}, // USE_WEFAC_IN_NETWORK + }, + }, }; return fully_supported_keywords_strings; @@ -66,6 +78,12 @@ const SupportedKeywordItems& fullySupported() { static const SupportedKeywordItems fully_supported_keywords_double = { + { + "NEFAC", + { + {2,{true, [](double x) { return x > 0 && x <= 1.0; }, "NEFAC(EFF_FACTOR: Efficiency must be in the range (0,1]"}}, // NETWORK_EFF_FACTOR + }, + }, { "WPIMULT", { diff --git a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp index 2dada5013..3e5fc99c0 100644 --- a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp @@ -111,12 +111,6 @@ partiallySupported() {8,{true, allow_values {"NO"}, "GECON(ENDRUN): End run not implemented"}}, }, }, - { - "GEFAC", - { - {3,{true, allow_values {"YES"}, "GEFAC(GRPNETWK): Extended Network Model efficiency NO option not implemented"}}, // TRANSFER_EXT_NET - }, - }, { "GRIDOPTS", { @@ -273,12 +267,6 @@ partiallySupported() {5,{true, allow_values {"NO"}, "WAGHYSTR(WATER_MODEL): only the NO option is supported – will STOP"}}, // WATER_MODEL }, }, - { - "WEFAC", - { - {3,{true, allow_values {"YES"}, "WEFAC(WELNETWK): only the YES option is supported"}}, // EXTENDED_NETWORK_OPT - }, - }, { "WELSPECS", { diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp index da61bbd8b..ce6196ce3 100644 --- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp +++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp @@ -397,7 +397,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords() {"NARROW", {true, std::nullopt}}, {"NCONSUMP", {true, std::nullopt}}, {"NCOMPS", {false, std::nullopt}}, - {"NEFAC", {true, std::nullopt}}, {"NETCOMPA", {true, std::nullopt}}, {"NEXT", {false, std::nullopt}}, {"NEXTSTPL", {true, std::nullopt}}, diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index a13874990..c5b665ff6 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -118,13 +118,13 @@ template class WellContributions; // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag - typedef Dune::FieldVector VectorBlockType; - typedef Dune::BlockVector BVector; + using VectorBlockType = Dune::FieldVector; + using BVector = Dune::BlockVector; - typedef BlackOilPolymerModule PolymerModule; - typedef BlackOilMICPModule MICPModule; + using PolymerModule = BlackOilPolymerModule; + using MICPModule = BlackOilMICPModule; - // For the conversion between the surface volume rate and resrevoir voidage rate + // For the conversion between the surface volume rate and reservoir voidage rate using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index ca69376de..d1c3b48f6 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1373,6 +1373,10 @@ updateAndCommunicateGroupData(const int reportStepIdx, reportStepIdx, well_state_nupcol, this->groupState()); + WellGroupHelpers::updateNetworkLeafNodeProductionRates(schedule(), + reportStepIdx, + well_state_nupcol, + this->groupState()); WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index f1a78796f..73643c69a 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -41,6 +41,7 @@ GroupState GroupState::serializationTestObject() { GroupState result(3); result.m_production_rates = {{"test1", {1.0, 2.0}}}; + result.m_network_leaf_node_production_rates={{"test1", {1.0, 20}}}; result.production_controls = {{"test2", Group::ProductionCMode::LRAT}}; result.prod_red_rates = {{"test3", {3.0, 4.0, 5.0}}}; result.inj_red_rates = {{"test4", {6.0, 7.0}}}; @@ -61,6 +62,7 @@ template bool GroupState::operator==(const GroupState& other) const { return this->m_production_rates == other.m_production_rates && + this->m_network_leaf_node_production_rates == other.m_network_leaf_node_production_rates && this->production_controls == other.production_controls && this->prod_red_rates == other.prod_red_rates && this->inj_red_rates == other.inj_red_rates && @@ -93,6 +95,16 @@ void GroupState::update_production_rates(const std::string& gname, this->m_production_rates[gname] = rates; } +template +void GroupState::update_network_leaf_node_production_rates(const std::string& gname, + const std::vector& rates) +{ + if (rates.size() != this->num_phases) + throw std::logic_error("Wrong number of phases"); + + this->m_network_leaf_node_production_rates[gname] = rates; +} + template const std::vector& GroupState::production_rates(const std::string& gname) const @@ -104,6 +116,17 @@ GroupState::production_rates(const std::string& gname) const return group_iter->second; } +template +const std::vector& +GroupState::network_leaf_node_production_rates(const std::string& gname) const +{ + auto group_iter = this->m_network_leaf_node_production_rates.find(gname); + if (group_iter == this->m_network_leaf_node_production_rates.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} + //------------------------------------------------------------------------- template diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 06bf1466f..41cc270bd 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -52,7 +52,10 @@ public: bool has_production_rates(const std::string& gname) const; void update_production_rates(const std::string& gname, const std::vector& rates); + void update_network_leaf_node_production_rates(const std::string& gname, + const std::vector& rates); const std::vector& production_rates(const std::string& gname) const; + const std::vector& network_leaf_node_production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); Scalar well_group_thp(const std::string& gname) const; @@ -130,6 +133,7 @@ public: auto forAllGroupData = [&](auto& func) { iterateContainer(m_production_rates, func); + iterateContainer(m_network_leaf_node_production_rates, func); iterateContainer(prod_red_rates, func); iterateContainer(inj_red_rates, func); iterateContainer(inj_resv_rates, func); @@ -187,6 +191,7 @@ public: { serializer(num_phases); serializer(m_production_rates); + serializer(m_network_leaf_node_production_rates); serializer(production_controls); serializer(group_thp); serializer(prod_red_rates); @@ -205,6 +210,7 @@ public: private: std::size_t num_phases{}; std::map> m_production_rates; + std::map> m_network_leaf_node_production_rates; std::map production_controls; std::map> prod_red_rates; std::map> inj_red_rates; diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 1c2a99507..233d9c48c 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -152,6 +152,7 @@ update(const WellState& well_state, } } } + init(); } template @@ -208,6 +209,7 @@ updateNewton(const BVectorWell& dwells, if (stop_or_zero_rate_target) { value_[0][WQTotal] = 0.; } + init(); } template @@ -515,7 +517,7 @@ template typename MultisegmentWellPrimaryVariables::EvalWell MultisegmentWellPrimaryVariables:: volumeFraction(const int seg, - const unsigned compIdx) const + const int compIdx) const { if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { return evaluation_[seg][WFrac]; @@ -584,7 +586,7 @@ typename MultisegmentWellPrimaryVariables::EvalWell MultisegmentWellPrimaryVariables:: getSegmentRateUpwinding(const int seg, const int seg_upwind, - const std::size_t comp_idx) const + const int comp_idx) const { // the result will contain the derivative with respect to WQTotal in segment seg, // and the derivatives with respect to WFrac GFrac in segment seg_upwind. diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 3dd05a61c..951d34b4c 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -115,7 +115,7 @@ public: //! \brief Returns upwinding rate for a component in a segment. EvalWell getSegmentRateUpwinding(const int seg, const int seg_upwind, - const std::size_t comp_idx) const; + const int comp_idx) const; //! \brief Get bottomhole pressure. EvalWell getBhp() const; @@ -154,7 +154,7 @@ private: //! \brief Returns volume fraction for component in a segment. EvalWell volumeFraction(const int seg, - const unsigned compIdx) const; + const int compIdx) const; //! \brief The values for the primary variables //! \details Based on different solution strategies, the wells can have different primary variables diff --git a/opm/simulators/wells/RatioCalculator.cpp b/opm/simulators/wells/RatioCalculator.cpp index d64990b7c..72d8ea5cd 100644 --- a/opm/simulators/wells/RatioCalculator.cpp +++ b/opm/simulators/wells/RatioCalculator.cpp @@ -51,9 +51,9 @@ namespace Opm { template RatioCalculator:: -RatioCalculator(unsigned gasCompIdx, - unsigned oilCompIdx, - unsigned waterCompIdx, +RatioCalculator(int gasCompIdx, + int oilCompIdx, + int waterCompIdx, std::string_view name) : gasComp_{gasCompIdx} , oilComp_(oilCompIdx) diff --git a/opm/simulators/wells/RatioCalculator.hpp b/opm/simulators/wells/RatioCalculator.hpp index 8121c3bb1..985881b88 100644 --- a/opm/simulators/wells/RatioCalculator.hpp +++ b/opm/simulators/wells/RatioCalculator.hpp @@ -39,9 +39,9 @@ class RatioCalculator public: using Scalar = decltype(getValue(Value{})); - RatioCalculator(unsigned gasCompIdx, - unsigned oilCompIdx, - unsigned waterCompIdx, + RatioCalculator(int gasCompIdx, + int oilCompIdx, + int waterCompIdx, std::string_view name); void disOilVapWatVolumeRatio(Value& volumeRatio, @@ -91,9 +91,9 @@ public: const bool isProducer) const; private: - unsigned gasComp_; - unsigned oilComp_; - unsigned waterComp_; + int gasComp_; + int oilComp_; + int waterComp_; std::string name_; }; diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 3b4c17f20..d333ee0fc 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -731,7 +731,7 @@ connectionRateFoam(const std::vector& cq_s, } case Phase::SOLVENT: { if constexpr (Indices::enableSolvent) - return static_cast(Indices::contiSolventEqIdx); + return Indices::contiSolventEqIdx; else OPM_DEFLOG_THROW(std::runtime_error, "Foam transport phase is SOLVENT but SOLVENT is not activated.", deferred_logger); } diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index 17410e2f3..b6b13bf01 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -223,6 +223,7 @@ update(const WellState& well_state, // BHP value_[Bhp] = ws.bhp; + init(); } template @@ -301,6 +302,7 @@ updateNewton(const BVectorWell& dwells, // so that bhp constaint can be an active control when needed. constexpr Scalar bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal; value_[Bhp] = std::max(value_[Bhp] - dx1_limited, bhp_lower_limit); + init(); } template @@ -320,6 +322,7 @@ updateNewtonPolyMW(const BVectorWell& dwells) value_[pskin_index] -= relaxation_factor * dx_pskin; } } + init(); } template @@ -446,7 +449,7 @@ copyToWellStatePolyMW(WellState& well_state) const template typename StandardWellPrimaryVariables::EvalWell StandardWellPrimaryVariables:: -volumeFraction(const unsigned compIdx) const +volumeFraction(const int compIdx) const { if (FluidSystem::numActivePhases() == 1) { return EvalWell(numWellEq_ + Indices::numEq, 1.0); @@ -456,7 +459,7 @@ volumeFraction(const unsigned compIdx) const return evaluation_[GFrac]; } - if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) { + if (Indices::enableSolvent && compIdx == Indices::contiSolventEqIdx) { return evaluation_[SFrac]; } diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index f482af12e..7c0c551ac 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -157,7 +157,7 @@ private: DeferredLogger& deferred_logger) const; //! \brief Returns volume fraction for a component. - EvalWell volumeFraction(const unsigned compIdx) const; + EvalWell volumeFraction(const int compIdx) const; //! \brief Handle non-reasonable fractions due to numerical overshoot. void processFractions(); diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 86e9ca791..74425033d 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -90,14 +90,15 @@ namespace Opm { const Opm::WellState& wellState, const int reportStepIdx, const int phasePos, - const bool injector) + const bool injector, + const bool network) { Scalar rate = 0.0; for (const std::string& groupName : group.groups()) { const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx); - const auto& gefac = groupTmp.getGroupEfficiencyFactor(); - rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector); + const auto& gefac = groupTmp.getGroupEfficiencyFactor(network); + rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector, network); } // only sum satelite production once @@ -126,8 +127,7 @@ namespace Opm { if (wellEcl.getStatus() == Opm::Well::Status::SHUT) continue; - const Scalar factor = wellEcl.getEfficiencyFactor() * - wellState[wellEcl.name()].efficiency_scaling_factor; + const Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor; const auto& ws = wellState.well(well_index.value()); if (res_rates) { const auto& well_rates = ws.reservoir_rates; @@ -744,6 +744,30 @@ updateGroupProductionRates(const Group& group, group_state.update_production_rates(group.name(), rates); } +template +void WellGroupHelpers:: +updateNetworkLeafNodeProductionRates(const Schedule& schedule, + const int reportStepIdx, + const WellState& wellState, + GroupState& group_state) +{ + const auto& network = schedule[reportStepIdx].network(); + if (network.active()) { + const int np = wellState.numPhases(); + for (const auto& group_name : network.leaf_nodes()) { + assert(schedule[reportStepIdx].groups.has(group_name)); + const auto& group = schedule[reportStepIdx].groups.get(group_name); + std::vector network_rates(np, 0.0); + if (group.numWells() > 0) { + for (int phase = 0; phase < np; ++phase) { + network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true); + } + } + group_state.update_network_leaf_node_production_rates(group_name, network_rates); + } + } +} + template void WellGroupHelpers:: updateREINForGroups(const Group& group, @@ -928,23 +952,15 @@ computeNetworkPressures(const Network::ExtNetwork& network, continue; } - node_inflows[node] = group_state.production_rates(node); + node_inflows[node] = group_state.network_leaf_node_production_rates(node); // Add the ALQ amounts to the gas rates if requested. if (network.node(node).add_gas_lift_gas()) { const auto& group = schedule.getGroup(node, report_time_step); for (const std::string& wellname : group.wells()) { const Well& well = schedule.getWell(wellname, report_time_step); + if (well.isInjector() || !well_state.isOpen(wellname)) continue; - // Here we use the efficiency unconditionally, but if WEFAC item 3 - // for the well is false (it defaults to true) then we should NOT use - // the efficiency factor. Fixing this requires not only changing the - // code here, but also: - // - Adding a member to the well for this flag, and setting it in Schedule::handleWEFAC(). - // - Making the wells' maximum flows (i.e. not time-averaged by using a efficiency factor) - // available and using those (for wells with WEFAC(3) true only) when accumulating group - // rates, but ONLY for network calculations. - const Scalar efficiency = well.getEfficiencyFactor() * - well_state.getGlobalEfficiencyScalingFactor(wellname); + const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) * well_state.getGlobalEfficiencyScalingFactor(wellname); node_inflows[node][BlackoilPhases::Vapour] += well_state.getALQ(wellname) * efficiency; } } @@ -962,13 +978,14 @@ computeNetworkPressures(const Network::ExtNetwork& network, // Add downbranch rates to upbranch. std::vector& up = node_inflows[(*upbranch).uptree_node()]; const std::vector& down = node_inflows[node]; + // We now also support NEFAC + const Scalar efficiency = network.node(node).efficiency(); if (up.empty()) { - up = down; - } else { - assert (up.size() == down.size()); - for (std::size_t ii = 0; ii < up.size(); ++ii) { - up[ii] += down[ii]; - } + up = std::vector(down.size(), 0.0); + } + assert (up.size() == down.size()); + for (std::size_t ii = 0; ii < up.size(); ++ii) { + up[ii] += efficiency*down[ii]; } } } diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 62bd15b28..c5d0e23ca 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -56,7 +56,8 @@ public: const Opm::WellState& wellState, const int reportStepIdx, const int phasePos, - const bool injector); + const bool injector, + const bool network = false); static Scalar satelliteProduction(const ScheduleState& sched, const std::vector& groups, @@ -199,6 +200,12 @@ public: const WellState& wellState, GroupState& group_state); + static void updateNetworkLeafNodeProductionRates(const Schedule& schedule, + const int reportStepIdx, + const WellState& wellState, + GroupState& group_state); + + static void updateWellRatesFromGroupTargetScale(const Scalar scale, const Group& group, const Schedule& schedule, diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index 1150bdc18..983318488 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -78,7 +78,7 @@ flowPhaseToModelCompIdx(const int phaseIdx) const template int WellInterfaceIndices:: -modelCompIdxToFlowCompIdx(const unsigned compIdx) const +modelCompIdxToFlowCompIdx(const int compIdx) const { const auto& pu = this->phaseUsage(); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx) diff --git a/opm/simulators/wells/WellInterfaceIndices.hpp b/opm/simulators/wells/WellInterfaceIndices.hpp index 6ac6138e1..16146c9a9 100644 --- a/opm/simulators/wells/WellInterfaceIndices.hpp +++ b/opm/simulators/wells/WellInterfaceIndices.hpp @@ -41,7 +41,7 @@ public: using ModelParameters = typename WellInterfaceFluidSystem::ModelParameters; int flowPhaseToModelCompIdx(const int phaseIdx) const; - int modelCompIdxToFlowCompIdx(const unsigned compIdx) const; + int modelCompIdxToFlowCompIdx(const int compIdx) const; Scalar scalingFactor(const int phaseIdx) const; template