diff --git a/CMakeLists.txt b/CMakeLists.txt index 976f72086d3..ba26576e377 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -814,6 +814,62 @@ macro(opm-simulators_targets_hook) OPM_COMPILE_COMPONENTS_TEMPLATE_LIST=${OPM_COMPILE_COMPONENTS_TEMPLATE_LIST} ) + opm_add_executable( + TARGET + flowexp_comp2_2p + SOURCES + flowexperimental/comp/flowexp_comp2_2p_main.cpp + $ + LIBRARIES + opmsimulators + ) + + opm_add_executable( + TARGET + flowexp_comp3_2p + SOURCES + flowexperimental/comp/flowexp_comp3_2p_main.cpp + $ + LIBRARIES + opmsimulators + ) + + opm_add_test(flow_comp + ONLY_COMPILE + ALWAYS_ENABLE + DEPENDS + opmsimulators + LIBRARIES + opmsimulators + SOURCES + flowexperimental/comp/flow_comp.cpp + $ + ) + target_compile_definitions(flow_comp + PRIVATE + OPM_COMPILE_COMPONENTS_TEMPLATE_LIST=${OPM_COMPILE_COMPONENTS_TEMPLATE_LIST} + ) + + opm_add_executable( + TARGET + flow_comp2_2p + SOURCES + flowexperimental/comp/flow_comp2_2p.cpp + $ + LIBRARIES + opmsimulators + ) + + opm_add_executable( + TARGET + flow_comp3_2p + SOURCES + flowexperimental/comp/flow_comp3_2p.cpp + $ + LIBRARIES + opmsimulators + ) + if(BUILD_FLOW_FLOAT_VARIANTS) opm_add_executable( TARGET diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index cc4b74daa90..a41f7f92965 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -73,6 +73,7 @@ list (APPEND MAIN_SOURCE_FILES flowexperimental/BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp flowexperimental/comp/EmptyModel.hpp + flowexperimental/comp/flow_comp.hpp flowexperimental/comp/flowexp_comp.hpp flowexperimental/comp/wells/CompWellModel.hpp flowexperimental/comp/wells/CompWellModel_impl.hpp @@ -156,12 +157,14 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/MixingRateControls.cpp opm/simulators/flow/NlddReporting.cpp opm/simulators/flow/NonlinearSolver.cpp + opm/simulators/flow/NonlinearSystemCompositional.hpp + opm/simulators/flow/NonlinearSystemCompositional_impl.hpp opm/simulators/flow/partitionCells.cpp opm/simulators/flow/RFTContainer.cpp opm/simulators/flow/RSTConv.cpp opm/simulators/flow/RegionPhasePVAverage.cpp opm/simulators/flow/SimulatorConvergenceOutput.cpp - opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp + opm/simulators/flow/SimulatorFullyImplicit.cpp opm/simulators/flow/SimulatorReportBanners.cpp opm/simulators/flow/SimulatorSerializer.cpp opm/simulators/flow/SolutionContainers.cpp @@ -203,6 +206,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/timestepping/TimeStepControl.cpp opm/simulators/timestepping/gatherConvergenceReport.cpp opm/simulators/utils/ComponentName.cpp + opm/simulators/utils/ComponentNameFlowExp.cpp opm/simulators/utils/DeferredLogger.cpp opm/simulators/utils/FullySupportedFlowKeywords.cpp opm/simulators/utils/ParallelFileMerger.cpp @@ -938,13 +942,17 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/Banners.hpp opm/simulators/flow/BaseAquiferModel.hpp opm/simulators/flow/BioeffectsContainer.hpp - opm/simulators/flow/BlackoilModel.hpp - opm/simulators/flow/BlackoilModel_impl.hpp + opm/simulators/flow/NonlinearSystemBlackOilReservoir.hpp + opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp opm/simulators/flow/BlackoilModelConvergenceMonitor.hpp - opm/simulators/flow/BlackoilModelNldd.hpp + opm/simulators/flow/NonlinearSystemNldd.hpp opm/simulators/flow/BlackoilModelParameters.hpp opm/simulators/flow/BlackoilModelProperties.hpp - opm/simulators/flow/BlackoilModelTPSA.hpp + opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp + opm/simulators/flow/NonlinearSystem.hpp + opm/simulators/flow/NonlinearSystem_impl.hpp + opm/simulators/flow/NonlinearSystemCompositional.hpp + opm/simulators/flow/NonlinearSystemCompositional_impl.hpp opm/simulators/flow/CO2H2Container.hpp opm/simulators/flow/CollectDataOnIORank.hpp opm/simulators/flow/CollectDataOnIORank_impl.hpp @@ -1013,8 +1021,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/RSTConv.hpp opm/simulators/flow/RegionPhasePVAverage.hpp opm/simulators/flow/SimulatorConvergenceOutput.hpp - opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp - opm/simulators/flow/SimulatorFullyImplicitBlackoil_impl.hpp + opm/simulators/flow/SimulatorFullyImplicit.hpp + opm/simulators/flow/SimulatorFullyImplicit_impl.hpp opm/simulators/flow/SimulatorReportBanners.hpp opm/simulators/flow/SimulatorSerializer.hpp opm/simulators/flow/SolutionContainers.hpp diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 0f624e104bd..1308c75db24 100644 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -81,6 +81,7 @@ function(add_test_compareECLFiles) FILENAME DEV_SIMULATOR SIMULATOR + REFERENCE_SIMULATOR ABS_TOL REL_TOL DIR @@ -112,7 +113,11 @@ function(add_test_compareECLFiles) if(PARAM_RESTART_SCHED STREQUAL "false" OR PARAM_RESTART_SCHED STREQUAL "true") list(APPEND DRIVER_ARGS -h ${PARAM_RESTART_SCHED}) endif() - list(APPEND DRIVER_ARGS -u ${PARAM_SIMULATOR}) + set(reference_simulator ${PARAM_SIMULATOR}) + if(PARAM_REFERENCE_SIMULATOR) + set(reference_simulator ${PARAM_REFERENCE_SIMULATOR}) + endif() + list(APPEND DRIVER_ARGS -u ${reference_simulator}) if(USE_DEV_SIMULATOR_IN_TESTS AND PARAM_DEV_SIMULATOR) set(PARAM_SIMULATOR ${PARAM_DEV_SIMULATOR}) endif() @@ -151,6 +156,7 @@ function(add_test_compareSeparateECLFiles) DIR2 DEV_SIMULATOR SIMULATOR + REFERENCE_SIMULATOR ABS_TOL REL_TOL IGNORE_EXTRA_KW @@ -187,6 +193,9 @@ function(add_test_compareSeparateECLFiles) if(PARAM_IGNORE_EXTRA_KW) list(APPEND DRIVER_ARGS -y ${PARAM_IGNORE_EXTRA_KW}) endif() + if(PARAM_REFERENCE_SIMULATOR) + list(APPEND DRIVER_ARGS -l ${PARAM_REFERENCE_SIMULATOR}) + endif() opm_add_test(${PARAM_PREFIX}_${PARAM_SIMULATOR}+${PARAM_CASENAME} EXE_TARGET ${PARAM_SIMULATOR} @@ -676,6 +685,16 @@ add_test_runSimulator( --enable-vtk-output=true ) +add_test_runSimulator(CASENAME 1dcompositional_flowexp_comp3_2p + FILENAME 1D_COMP + SIMULATOR flowexp_comp3_2p + DIR compositional) + +add_test_runSimulator(CASENAME 1dcompositional_flow_comp3_2p + FILENAME 1D_COMP + SIMULATOR flow_comp3_2p + DIR compositional) + # Tests that are run based on simulator results, but not necessarily direct comparison to reference results add_test_runSimulator( CASENAME diff --git a/comparisonTests.cmake b/comparisonTests.cmake index 1e95e42d272..57068355843 100644 --- a/comparisonTests.cmake +++ b/comparisonTests.cmake @@ -31,3 +31,29 @@ add_test_compareSeparateECLFiles( MPI_PROCS 1 ) + +add_test_compareSeparateECLFiles( + CASENAME + 1dcompositional_flow_comp3_2p_vs_flowexp_comp3_2p + DIR1 + compositional + FILENAME1 + 1D_COMP + DIR2 + compositional + FILENAME2 + 1D_COMP + SIMULATOR + flow_comp3_2p + REFERENCE_SIMULATOR + flowexp_comp3_2p + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + TEST_ARGS + --solver-max-time-step-in-days=0.4 + MPI_PROCS + 1 +) + diff --git a/flow/flow_biofilm.cpp b/flow/flow_biofilm.cpp index 4b650f98199..a84138bd789 100644 --- a/flow/flow_biofilm.cpp +++ b/flow/flow_biofilm.cpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_blackoil.cpp b/flow/flow_blackoil.cpp index ce284beb32b..78b2606bbd5 100644 --- a/flow/flow_blackoil.cpp +++ b/flow/flow_blackoil.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_blackoil_legacyassembly.cpp b/flow/flow_blackoil_legacyassembly.cpp index 00ddd9e912b..fdf45f71cb3 100644 --- a/flow/flow_blackoil_legacyassembly.cpp +++ b/flow/flow_blackoil_legacyassembly.cpp @@ -19,7 +19,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_blackoil_nohyst.cpp b/flow/flow_blackoil_nohyst.cpp index 457cbaeaeb0..2de2fccc769 100644 --- a/flow/flow_blackoil_nohyst.cpp +++ b/flow/flow_blackoil_nohyst.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_blackoil_temp.cpp b/flow/flow_blackoil_temp.cpp index d02b56a4bbc..e5a93d9c58c 100644 --- a/flow/flow_blackoil_temp.cpp +++ b/flow/flow_blackoil_temp.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_blackoil_tpsa.cpp b/flow/flow_blackoil_tpsa.cpp index 9f3a2015be3..e01bb622e33 100644 --- a/flow/flow_blackoil_tpsa.cpp +++ b/flow/flow_blackoil_tpsa.cpp @@ -73,7 +73,7 @@ struct Problem template struct NonlinearSystem -{ using type = BlackoilModelTPSA; }; +{ using type = NonlinearSystemBlackOilReservoirTPSA; }; } // namespace Opm::Properties diff --git a/flow/flow_brine.cpp b/flow/flow_brine.cpp index b1e7a029a32..b8dcc4bcec6 100644 --- a/flow/flow_brine.cpp +++ b/flow/flow_brine.cpp @@ -23,7 +23,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_brine_precsalt_vapwat.cpp b/flow/flow_brine_precsalt_vapwat.cpp index e03d3da6328..2152fa6b0c6 100644 --- a/flow/flow_brine_precsalt_vapwat.cpp +++ b/flow/flow_brine_precsalt_vapwat.cpp @@ -23,7 +23,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_brine_saltprecipitation.cpp b/flow/flow_brine_saltprecipitation.cpp index c1e3a1c5ce3..72e3db3a1fe 100644 --- a/flow/flow_brine_saltprecipitation.cpp +++ b/flow/flow_brine_saltprecipitation.cpp @@ -23,7 +23,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_energy.cpp b/flow/flow_energy.cpp index 650f024f53b..d140d155f21 100644 --- a/flow/flow_energy.cpp +++ b/flow/flow_energy.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include #include #include diff --git a/flow/flow_extbo.cpp b/flow/flow_extbo.cpp index 7c9340240ce..cb03d1fa8ee 100644 --- a/flow/flow_extbo.cpp +++ b/flow/flow_extbo.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_foam.cpp b/flow/flow_foam.cpp index 8e1109a165b..0f17c2904e4 100644 --- a/flow/flow_foam.cpp +++ b/flow/flow_foam.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_gasoil.cpp b/flow/flow_gasoil.cpp index 1ab5e63a241..2f9f2f3bed9 100644 --- a/flow/flow_gasoil.cpp +++ b/flow/flow_gasoil.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gasoil_energy.cpp b/flow/flow_gasoil_energy.cpp index 0a826d54132..0dedc161b02 100644 --- a/flow/flow_gasoil_energy.cpp +++ b/flow/flow_gasoil_energy.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include #include diff --git a/flow/flow_gasoildiffuse.cpp b/flow/flow_gasoildiffuse.cpp index 1177d18ac42..6b0f6cb8459 100644 --- a/flow/flow_gasoildiffuse.cpp +++ b/flow/flow_gasoildiffuse.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include #include diff --git a/flow/flow_gaswater.cpp b/flow/flow_gaswater.cpp index 0d9aec6d9d8..010b3703bd9 100644 --- a/flow/flow_gaswater.cpp +++ b/flow/flow_gaswater.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_brine.cpp b/flow/flow_gaswater_brine.cpp index 1062f4c95fd..752788463b7 100644 --- a/flow/flow_gaswater_brine.cpp +++ b/flow/flow_gaswater_brine.cpp @@ -24,7 +24,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_dissolution.cpp b/flow/flow_gaswater_dissolution.cpp index 393b89f441d..a258bf9a380 100644 --- a/flow/flow_gaswater_dissolution.cpp +++ b/flow/flow_gaswater_dissolution.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_dissolution_diffuse.cpp b/flow/flow_gaswater_dissolution_diffuse.cpp index 930519fa425..e50233e3f45 100644 --- a/flow/flow_gaswater_dissolution_diffuse.cpp +++ b/flow/flow_gaswater_dissolution_diffuse.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_dissolution_tpsa.cpp b/flow/flow_gaswater_dissolution_tpsa.cpp index c9d8094a26e..83870bd4083 100644 --- a/flow/flow_gaswater_dissolution_tpsa.cpp +++ b/flow/flow_gaswater_dissolution_tpsa.cpp @@ -103,7 +103,7 @@ struct Problem template struct NonlinearSystem -{ using type = BlackoilModelTPSA; }; +{ using type = NonlinearSystemBlackOilReservoirTPSA; }; } // namespace Opm::Properties diff --git a/flow/flow_gaswater_energy.cpp b/flow/flow_gaswater_energy.cpp index a35dda24aae..25b9f1d1026 100644 --- a/flow/flow_gaswater_energy.cpp +++ b/flow/flow_gaswater_energy.cpp @@ -21,7 +21,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_saltprec_energy.cpp b/flow/flow_gaswater_saltprec_energy.cpp index 7a712deac69..2baadda48bf 100644 --- a/flow/flow_gaswater_saltprec_energy.cpp +++ b/flow/flow_gaswater_saltprec_energy.cpp @@ -24,7 +24,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_saltprec_vapwat.cpp b/flow/flow_gaswater_saltprec_vapwat.cpp index c69b7cbd5a1..fb45830fa06 100644 --- a/flow/flow_gaswater_saltprec_vapwat.cpp +++ b/flow/flow_gaswater_saltprec_vapwat.cpp @@ -24,7 +24,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_gaswater_solvent.cpp b/flow/flow_gaswater_solvent.cpp index 1ee649ac7cd..6a6f5e75984 100644 --- a/flow/flow_gaswater_solvent.cpp +++ b/flow/flow_gaswater_solvent.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_micp.cpp b/flow/flow_micp.cpp index defb243ee49..b7fd1b592e3 100644 --- a/flow/flow_micp.cpp +++ b/flow/flow_micp.cpp @@ -26,7 +26,7 @@ #include #include -#include +#include namespace Opm { /*! diff --git a/flow/flow_oilwater.cpp b/flow/flow_oilwater.cpp index 6a7ef616ee1..64084531de3 100644 --- a/flow/flow_oilwater.cpp +++ b/flow/flow_oilwater.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include #include diff --git a/flow/flow_oilwater_brine.cpp b/flow/flow_oilwater_brine.cpp index b997edb92fe..d08b67211d8 100644 --- a/flow/flow_oilwater_brine.cpp +++ b/flow/flow_oilwater_brine.cpp @@ -24,7 +24,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_oilwater_polymer.cpp b/flow/flow_oilwater_polymer.cpp index 34a78016668..e800d7404fd 100644 --- a/flow/flow_oilwater_polymer.cpp +++ b/flow/flow_oilwater_polymer.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_oilwater_polymer_injectivity.cpp b/flow/flow_oilwater_polymer_injectivity.cpp index 41ac7294287..4b1512fea8d 100644 --- a/flow/flow_oilwater_polymer_injectivity.cpp +++ b/flow/flow_oilwater_polymer_injectivity.cpp @@ -22,7 +22,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_onephase_tpsa.cpp b/flow/flow_onephase_tpsa.cpp index deae40ba7a2..e848fd09c43 100644 --- a/flow/flow_onephase_tpsa.cpp +++ b/flow/flow_onephase_tpsa.cpp @@ -94,7 +94,7 @@ struct Problem template struct NonlinearSystem -{ using type = BlackoilModelTPSA; }; +{ using type = NonlinearSystemBlackOilReservoirTPSA; }; } // namespace Opm::Properties diff --git a/flow/flow_polymer.cpp b/flow/flow_polymer.cpp index 3456b41321c..4b81a8ed579 100644 --- a/flow/flow_polymer.cpp +++ b/flow/flow_polymer.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_solvent.cpp b/flow/flow_solvent.cpp index f331c854469..4ed590038a9 100644 --- a/flow/flow_solvent.cpp +++ b/flow/flow_solvent.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flow/flow_solvent_foam.cpp b/flow/flow_solvent_foam.cpp index eadd1b740ec..887e726b485 100644 --- a/flow/flow_solvent_foam.cpp +++ b/flow/flow_solvent_foam.cpp @@ -20,7 +20,7 @@ #include #include -#include +#include #include namespace Opm { diff --git a/flowexperimental/comp/EmptyModel.hpp b/flowexperimental/comp/EmptyModel.hpp index 9f0c75efcc2..81b000bfdd4 100644 --- a/flowexperimental/comp/EmptyModel.hpp +++ b/flowexperimental/comp/EmptyModel.hpp @@ -58,6 +58,10 @@ class EmptyModel : public BaseAuxiliaryModule template void deserialize(Restarter& /*res*/){}; + template + void serializeOp(Serializer&) + {} + void beginEpisode(){}; void beginTimeStep(){}; void beginIteration(){}; diff --git a/flowexperimental/comp/flow_comp.cpp b/flowexperimental/comp/flow_comp.cpp new file mode 100644 index 00000000000..1a334da4468 --- /dev/null +++ b/flowexperimental/comp/flow_comp.cpp @@ -0,0 +1,119 @@ +/* + 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 + +#include +#include + +#include + +#include + +#include +#include +#include + +#include "flow_comp.hpp" + +template +std::tuple +runComponent(int runtimeComponent, bool water, int argc, char** argv) +{ + if (runtimeComponent == compileTimeComponent) { + if (water) + return std::make_tuple(true, Opm::dispatchFlowComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowComp(argc, argv)); + } + return std::make_tuple(false, EXIT_FAILURE); +} + + +/** + * @brief Runs a specified runtime component. + * + * This function checks if the provided runtime component matches the current compile-time component. + * If they match, it dispatches the simulator for the number of components and returns + * a tuple where the second element is the result of the execution. + * + * Otherwise, it recursively calls itself with the next component in the list. + * + * @param runtimecomponent The runtime component identifier to be executed. + * @param argc The number of command-line arguments. + * @param argv The array of command-line arguments. + * @return A tuple containing a boolean indicating if the component was executed and an integer result of the execution. + * + * @note We have two non-variadic templates to be able to properly overload for the base case. + */ +template +std::tuple +runComponent(int runtimecomponent, bool water, int argc, char** argv) +{ + if (currentCompileTimeComponent == runtimecomponent) { + if (water) + return std::make_tuple(true, Opm::dispatchFlowComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowComp(argc, argv)); + } + return runComponent(runtimecomponent, water, argc, argv); +} + +int +main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::FlowCompProblem<0, true>; + Opm::registerEclTimeSteppingParameters(); + + // At the moment, this is probably as optimal as can be. + // We only read the RUNSPEC of the Deck file to get the numComp, + // and for this we need to first read the CLI arguments. + Opm::setupParameters_(argc, const_cast(argv), true, false, true, 0); + + auto inputFilename + = Opm::FlowGenericVanguard::canonicalDeckPath(Opm::Parameters::Get()); + + // Only read the RUNSPEC section of the deck + const auto deck + = 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); + + Opm::Parameters::reset(); + + auto [componentSupported, executionStatus] + = runComponent(numComps, wat, argc, argv); + + if (!componentSupported) { + fmt::print("Deck has {} components, not supported. In this build of the simulator, we support the " + "following number of components:\n\n\t{}\n\nNote that the supported components can be changed " + "when configuring CMake through the OPM_COMPILE_COMPONENTS option. That is, run cmake as " + "eg\n\n\tcmake ... -DOPM_COMPILE_COMPONENTS=\"4;7\"\n\nto compile only components 4 and 7, " + "or\n\n\tcmake ... -DOPM_COMPILE_COMPONENTS=3\n\nto compile only component 3, then " + "recompile.\n\nExiting.\n", + numComps, + fmt::join(std::array {OPM_COMPILE_COMPONENTS_TEMPLATE_LIST}, ", ")); + } + return executionStatus; +} diff --git a/flowexperimental/comp/flow_comp.hpp b/flowexperimental/comp/flow_comp.hpp new file mode 100644 index 00000000000..656f7ecfb39 --- /dev/null +++ b/flowexperimental/comp/flow_comp.hpp @@ -0,0 +1,296 @@ +/* + 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 . +*/ +#ifndef FLOW_COMP_HPP +#define FLOW_COMP_HPP + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +// // the current code use eclnewtonmethod adding other conditions to proceed_ should do the trick for KA +// // adding linearshe sould be chaning the update_ function in the same class with condition that the error is reduced. +// the trick is to be able to recalculate the residual from here. +// unsure where the timestepping is done from suggestedtime?? +// suggestTimeStep is taken from newton solver in problem.limitTimestep +namespace Opm { + +template +int dispatchFlowComp(int argc, char** argv); + +} // namespace Opm + +namespace Opm::Properties { + +template +struct DiscNewtonMethod; + +template +struct WellModel; + +template +struct TracerModel; + +namespace TTag { + +template +struct FlowCompProblem { + using InheritsFrom = std::tuple; +}; + +} + +template +struct SparseMatrixAdapter> +{ +private: + using Scalar = GetPropType; + enum { numEq = getPropValue() }; + using Block = MatrixBlock; + +public: + using type = typename Linear::IstlSparseMatrixAdapter; +}; + +// NOTE: Do NOT override NewtonMethod here. +// The inherited model configuration decides the solver implementation. + +#if 0 +template +struct SolidEnergyLaw +{ +private: + using Scalar = GetPropType; + using FluidSystem = GetPropType; + +public: + using EclThermalLawManager = ::Opm::EclThermalLawManager; + + using type = typename EclThermalLawManager::SolidEnergyLaw; +}; + +// Set the material law for thermal conduction +template +struct ThermalConductionLaw +{ +private: + using Scalar = GetPropType; + using FluidSystem = GetPropType; + +public: + using EclThermalLawManager = ::Opm::EclThermalLawManager; + + using type = typename EclThermalLawManager::ThermalConductionLaw; +}; + + +template +struct SpatialDiscretizationSplice +{ + using type = TTag::EcfvDiscretization; +}; + +template +struct LocalLinearizerSplice +{ + using type = TTag::AutoDiffLocalLinearizer; +}; +#endif + +// Set the problem property +template +struct Problem> +{ + using type = FlowProblemComp; +}; + +template +struct AquiferModel> { + using type = EmptyModel; +}; + +template +struct WellModel> { + using type = CompWellModel; +}; + +template +struct TracerModel> { + using type = EmptyModel; +}; + + +template +struct FlashSolver> { +private: + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Evaluation = GetPropType; + +public: + using type = Opm::PTFlash; +}; + + +template +struct NumComp { using type = UndefinedProperty; }; + +// TODO: this is unfortunate, have to check why we need to hard-code it +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; }; + + template + struct Temperature { + using type = GetPropType; + static constexpr type value = 423.25; + }; +#endif + +template +struct FluidSystem> +{ +private: + using Scalar = GetPropType; + static constexpr int num_comp = getPropValue(); + static constexpr bool enable_water = getPropValue(); + +public: + using type = Opm::GenericOilGasWaterFluidSystem; +}; +template +struct EnableMech> { + static constexpr bool value = false; +}; + +template +struct EnableDisgasInWater> { + static constexpr bool value = false; +}; + +template +struct Stencil> +{ +private: + using Scalar = GetPropType; + using GridView = GetPropType; + +public: + using type = EcfvStencil; +}; + +template +struct EnableApiTracking> { + static constexpr bool value = false; +}; + +template +struct EnableSaltPrecipitation> { + static constexpr bool value = false; +}; +template +struct EnablePolymerMW> { + static constexpr bool value = false; +}; + +template +struct EnablePolymer> { + static constexpr bool value = false; +}; + +template +struct EnableDispersion> { + static constexpr bool value = false; +}; + +template +struct EnableBrine> { + static constexpr bool value = false; +}; +template +struct EnableVapwat> { + static constexpr bool value = false; +}; + +template +struct EnableSolvent> { + static constexpr bool value = false; +}; +template +struct EnableFoam> { + static constexpr bool value = false; +}; +template +struct EnableExtbo> { + static constexpr bool value = false; +}; +template +struct EnableBioeffects> { + static constexpr bool value = false; +}; + +// disable thermal flux boundaries by default +#if 0 +template +struct EnableThermalFluxBoundaries { + static constexpr bool value = false; +}; +#endif + +} // namespace Opm::Properties + +namespace Opm { + +template +int dispatchFlowComp(int argc, char** argv) +{ + using TypeTag = Properties::TTag::FlowCompProblem; + auto mainObject = std::make_unique(argc, argv); + const auto ret = mainObject->runStatic(); + mainObject.reset(); + return ret; +} + +} // namespace Opm + +#endif diff --git a/flowexperimental/comp/flow_comp2_2p.cpp b/flowexperimental/comp/flow_comp2_2p.cpp new file mode 100644 index 00000000000..a819b1a54bb --- /dev/null +++ b/flowexperimental/comp/flow_comp2_2p.cpp @@ -0,0 +1,29 @@ +/* + 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 "flow_comp.hpp" + +int main(int argc, char** argv) +{ + return Opm::dispatchFlowComp<2, false>(argc, argv); +} diff --git a/flowexperimental/comp/flow_comp3_2p.cpp b/flowexperimental/comp/flow_comp3_2p.cpp new file mode 100644 index 00000000000..3c240636db2 --- /dev/null +++ b/flowexperimental/comp/flow_comp3_2p.cpp @@ -0,0 +1,29 @@ +/* + 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 "flow_comp.hpp" + +int main(int argc, char** argv) +{ + return Opm::dispatchFlowComp<3, false>(argc, argv); +} diff --git a/flowexperimental/comp/flowexp_comp2_2p_main.cpp b/flowexperimental/comp/flowexp_comp2_2p_main.cpp new file mode 100644 index 00000000000..7a9be82c398 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp2_2p_main.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" + +int main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<2, false>; + + Opm::registerEclTimeSteppingParameters(); + Opm::setupParameters_(argc, const_cast(argv), true, false, true, 0); + + return Opm::start(argc, argv, false); +} diff --git a/flowexperimental/comp/flowexp_comp3_2p_main.cpp b/flowexperimental/comp/flowexp_comp3_2p_main.cpp new file mode 100644 index 00000000000..49dc4101fb0 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp3_2p_main.cpp @@ -0,0 +1,39 @@ +/* + 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 + +#include "flowexp_comp.hpp" + +int main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<3, false>; + + Opm::Parameters::Register + ("Use adaptive time stepping between report steps"); + Opm::registerEclTimeSteppingParameters(); + Opm::setupParameters_(argc, const_cast(argv), true, false, true, 0); + + return Opm::start(argc, argv, false); +} diff --git a/flowexperimental/comp/wells/CompWellModel.hpp b/flowexperimental/comp/wells/CompWellModel.hpp index 5bf9d82be07..f657d172c1c 100644 --- a/flowexperimental/comp/wells/CompWellModel.hpp +++ b/flowexperimental/comp/wells/CompWellModel.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -35,7 +36,10 @@ #include +#include + +#include #include namespace Opm { @@ -44,7 +48,7 @@ class Schedule; struct NewtonIterationContext; template -class CompWellModel : WellConnectionAuxiliaryModule> +class CompWellModel : public WellConnectionAuxiliaryModule> { using ElementContext = GetPropType; using FluidSystem = GetPropType; @@ -91,15 +95,20 @@ class CompWellModel : WellConnectionAuxiliaryModule + void serializeOp(Serializer&) + {} + void beginEpisode() { beginReportStep(simulator_.episodeIndex()); } void beginReportStep(unsigned report_step); void beginTimeStep(); void beginIteration(); + void restoreLastValidState(); void init(); void endIteration() const {} - void endTimeStep() {} + void endTimeStep(); void endEpisode() {} void computeTotalRatesForDof(RateVector& /*rate*/, unsigned /*globalIdx*/) const; @@ -126,6 +135,19 @@ class CompWellModel : WellConnectionAuxiliaryModule& wellOpenTimes() const { return well_open_times_; } + const std::map& wellCloseTimes() const { return well_close_times_; } + const WellGroupEvents& reportStepStartEvents() const { return report_step_start_events_; } + bool forceShutWellByName(const std::string& well_name, double simulation_time, bool dont_shut_grup_wells); + + template + void setReservoirCouplingSlave(ReservoirCouplingSlave*) {} + + template + void setReservoirCouplingMaster(ReservoirCouplingMaster*) {} + bool getWellConvergence() const; // the following functions are not used while added to avoid modifying WellConnectionAuxiliaryModule.hpp @@ -146,6 +168,9 @@ class CompWellModel : WellConnectionAuxiliaryModule comp_well_states_; + // saved state at beginning of report step, used to restore on failed timestep + CompWellState last_valid_comp_well_states_; + // this is needed for parallel running, not all the wells will be in the same process std::vector wells_ecl_; std::vector > well_connection_data_; @@ -155,6 +180,10 @@ class CompWellModel : WellConnectionAuxiliaryModule well_open_times_; + std::map well_close_times_; + WellGroupEvents report_step_start_events_{}; void createWellContainer(); void initWellContainer(); diff --git a/flowexperimental/comp/wells/CompWellModel_impl.hpp b/flowexperimental/comp/wells/CompWellModel_impl.hpp index 6090e5121b3..3299c379691 100644 --- a/flowexperimental/comp/wells/CompWellModel_impl.hpp +++ b/flowexperimental/comp/wells/CompWellModel_impl.hpp @@ -26,6 +26,8 @@ #include #endif +#include + #include #include #include @@ -42,6 +44,7 @@ CompWellModel::CompWellModel(Simulator& simulator, const NewtonIteratio , comm_(simulator.gridView().comm()) , comp_config_(ecl_state_.compositionalConfig()) , comp_well_states_(comp_config_) + , last_valid_comp_well_states_(comp_config_) { local_num_cells_ = simulator.gridView().size(0); } @@ -52,9 +55,35 @@ CompWellModel:: beginReportStep(unsigned report_step) { // TODO: not considering the parallel running yet + report_step_start_events_ = schedule_[report_step].wellgroup_events(); wells_ecl_ = schedule_.getWells(report_step); + + constexpr auto events_mask = ScheduleEvents::WELL_STATUS_CHANGE | + ScheduleEvents::REQUEST_OPEN_WELL | + ScheduleEvents::REQUEST_SHUT_WELL; + for (const auto& well_ecl : wells_ecl_) { + if (!well_ecl.hasConnections()) { + continue; + } + + if (!report_step_start_events_.hasEvent(well_ecl.name(), events_mask)) { + continue; + } + + if (well_ecl.getStatus() == WellStatus::OPEN) { + well_open_times_.insert_or_assign(well_ecl.name(), simulator_.time()); + well_close_times_.erase(well_ecl.name()); + } + else if (well_ecl.getStatus() == WellStatus::SHUT) { + well_close_times_.insert_or_assign(well_ecl.name(), simulator_.time()); + well_open_times_.erase(well_ecl.name()); + } + } + initWellConnectionData(); initWellState(); + // Save the initial accepted state for this report step. + last_valid_comp_well_states_.copyDynamicStateFrom(comp_well_states_); } template @@ -66,6 +95,24 @@ beginTimeStep() initWellContainer(); } +template +void +CompWellModel:: +restoreLastValidState() +{ + comp_well_states_.copyDynamicStateFrom(last_valid_comp_well_states_); +} + +template +void +CompWellModel:: +endTimeStep() +{ + // Persist the accepted well state so failed retries restart from the last + // successful timestep rather than from the beginning of the report step. + last_valid_comp_well_states_.copyDynamicStateFrom(comp_well_states_); +} + template void CompWellModel:: @@ -83,6 +130,12 @@ createWellContainer() const auto nw = wells_ecl_.size(); // not considering the parallel running yet well_container_.clear(); for (auto w = 0 * nw; w < nw; ++w) { + const auto& well_name = wells_ecl_[w].name(); + if (comp_well_states_.has(well_name) + && comp_well_states_[well_name].status == WellStatus::SHUT) { + continue; + } + well_container_.emplace_back(std::make_shared>(wells_ecl_[w], w, well_connection_data_[w])); } } @@ -180,11 +233,13 @@ initWellState() OPM_END_PARALLEL_TRY_CATCH("ComposotionalWellModel::initializeWellState() failed: ", this->simulator_.vanguard().grid().comm()); - /* TODO: no prev well state for now */ - // \Note:: we are not supporting dynamic temperature, so the temperature is constant, so we just pass in a single value + // Start each report step from freshly initialized schedule state. Retry + // recovery still comes from last_valid_comp_well_states_ via + // restoreLastValidState()/endTimeStep(). this->comp_well_states_.init(this->wells_ecl_, cell_pressure, cell_temperature[0], cell_mole_fractions, this->well_connection_data_, - this->summary_state_); + this->summary_state_, + nullptr);//better but change result &this->last_valid_comp_well_states_); } @@ -293,6 +348,38 @@ recoverWellSolutionAndUpdateWellState(const BVector& x) } } +template +bool +CompWellModel:: +forceShutWellByName(const std::string& well_name, + double simulation_time, + bool) +{ + int well_was_shut = 0; + + if (comp_well_states_.has(well_name)) { + auto& well_state = comp_well_states_[well_name]; + if (well_state.status != WellStatus::SHUT) { + well_state.status = WellStatus::SHUT; + well_close_times_.insert_or_assign(well_name, simulation_time); + well_open_times_.erase(well_name); + + if (last_valid_comp_well_states_.has(well_name)) { + last_valid_comp_well_states_[well_name].status = WellStatus::SHUT; + } + + std::erase_if(well_container_, + [&well_name](const auto& well) + { return well->name() == well_name; }); + + well_was_shut = 1; + } + } + + well_was_shut = comm_.max(well_was_shut); + return well_was_shut == 1; +} + template bool CompWellModel:: diff --git a/flowexperimental/comp/wells/CompWellState.hpp b/flowexperimental/comp/wells/CompWellState.hpp index bcd705ef8a6..4ac55539a89 100644 --- a/flowexperimental/comp/wells/CompWellState.hpp +++ b/flowexperimental/comp/wells/CompWellState.hpp @@ -57,8 +57,17 @@ class CompWellState SingleWellState& operator[](const std::string& well_name); + bool has(const std::string& well_name) const + { return wells_.has(well_name); } + data::Wells report() const; + // Copy dynamic well variables while keeping each instance's own config reference. + void copyDynamicStateFrom(const CompWellState& other) + { + wells_ = other.wells_; + } + private: WellContainer wells_; diff --git a/flowexperimental/comp/wells/CompWellState_impl.hpp b/flowexperimental/comp/wells/CompWellState_impl.hpp index c2231c108e8..ac1082ea8ae 100644 --- a/flowexperimental/comp/wells/CompWellState_impl.hpp +++ b/flowexperimental/comp/wells/CompWellState_impl.hpp @@ -35,12 +35,19 @@ init(const std::vector& wells_ecl, const std::vector>& cell_mole_fractions, const std::vector >& well_connection_data, const SummaryState& summary_state, - const CompWellState* /*prev_well_state*/) + const CompWellState* prev_well_state) { this->base_init(wells_ecl, cell_pressures, temperature, cell_mole_fractions, well_connection_data, summary_state); - // TODO: for the simple case we have, I think we can just copy from the prev_well_state - // let us see how we gonna use it though + if (!prev_well_state) { + return; + } + + for (const auto& well_name : this->wells_.wells()) { + if (prev_well_state->wells_.has(well_name)) { + this->wells_[well_name].copyRuntimeStateFrom(prev_well_state->wells_[well_name]); + } + } } template diff --git a/flowexperimental/comp/wells/CompWell_impl.hpp b/flowexperimental/comp/wells/CompWell_impl.hpp index a699d2d5c9e..a2e96103148 100644 --- a/flowexperimental/comp/wells/CompWell_impl.hpp +++ b/flowexperimental/comp/wells/CompWell_impl.hpp @@ -17,6 +17,9 @@ along with OPM. If not, see . */ +#include +#include + #include #include @@ -563,6 +566,8 @@ updateWellControl(const SummaryState& summary_state, if (injection_controls.hasControl(Well::InjectorCMode::BHP) && current_control != WellInjectorCMode::BHP) { const Scalar bhp_limit = injection_controls.bhp_limit; const Scalar current_bhp = well_state.bhp; + OpmLog::debug(fmt::format("Well {} BHP control check: current_bhp={:.6e}, bhp_limit={:.6e}, exceeds_limit={}", + this->well_ecl_.name(), current_bhp, bhp_limit, (current_bhp > bhp_limit))); if (current_bhp > bhp_limit) { well_state.bhp = bhp_limit; well_state.injection_cmode = WellInjectorCMode::BHP; @@ -575,7 +580,12 @@ updateWellControl(const SummaryState& summary_state, // TODO: hack to get the injection rate const Scalar current_rate = std::accumulate(well_state.surface_phase_rates.begin(), well_state.surface_phase_rates.end(), 0.0); + OpmLog::debug(fmt::format("Well {} RATE control check: current_rate={:.6e}, rate_limit={:.6e}, bhp={:.6e}, phase_rates=[{}]", + this->well_ecl_.name(), current_rate, rate_limit, well_state.bhp, + fmt::join(well_state.surface_phase_rates, ", "))); if (current_rate > rate_limit) { + OpmLog::info(fmt::format("Well {} RATE control TRIGGERED: current_rate={:.6e} > rate_limit={:.6e}", + this->well_ecl_.name(), current_rate, rate_limit)); well_state.injection_cmode = WellInjectorCMode::RATE; changed = true; } diff --git a/flowexperimental/comp/wells/SingleCompWellState.hpp b/flowexperimental/comp/wells/SingleCompWellState.hpp index 093fceb84f8..e802bc452f7 100644 --- a/flowexperimental/comp/wells/SingleCompWellState.hpp +++ b/flowexperimental/comp/wells/SingleCompWellState.hpp @@ -94,6 +94,8 @@ class SingleCompWellState void update_injector_targets(const Well& well, const SummaryState& st); + void copyRuntimeStateFrom(const SingleCompWellState& other); + Scalar get_total_surface_rate() const; }; diff --git a/flowexperimental/comp/wells/SingleCompWellState_impl.hpp b/flowexperimental/comp/wells/SingleCompWellState_impl.hpp index acf8bcba0d7..2599f9c7ac8 100644 --- a/flowexperimental/comp/wells/SingleCompWellState_impl.hpp +++ b/flowexperimental/comp/wells/SingleCompWellState_impl.hpp @@ -145,6 +145,22 @@ update_producer_targets(const Well& well, } } +template +void SingleCompWellState:: +copyRuntimeStateFrom(const SingleCompWellState& other) +{ + // Keep the freshly initialized schedule-derived status, controls, and + // injector targets from base_init(); only reuse the dynamic state. + bhp = other.bhp; + surface_phase_rates = other.surface_phase_rates; + phase_fractions = other.phase_fractions; + reservoir_phase_rates = other.reservoir_phase_rates; + if (producer) { + total_molar_fractions = other.total_molar_fractions; + phase_molar_fractions = other.phase_molar_fractions; + } +} + template typename SingleCompWellState::Scalar SingleCompWellState:: diff --git a/opm/models/flash/flashprimaryvariables.hh b/opm/models/flash/flashprimaryvariables.hh index c61fcc02c86..e7d01c3a317 100644 --- a/opm/models/flash/flashprimaryvariables.hh +++ b/opm/models/flash/flashprimaryvariables.hh @@ -85,6 +85,14 @@ public: using ParentType::operator=; //!< Import base class assignment operators. + template + void serializeOp(Serializer& serializer) + { + for (unsigned i = 0; i < this->size(); ++i) { + serializer((*this)[i]); + } + } + /*! * \copydoc ImmisciblePrimaryVariables::assignMassConservative */ diff --git a/opm/models/nonlinear/newtonmethod.hh b/opm/models/nonlinear/newtonmethod.hh index 40a08987a80..abb4a25481f 100644 --- a/opm/models/nonlinear/newtonmethod.hh +++ b/opm/models/nonlinear/newtonmethod.hh @@ -192,6 +192,21 @@ public: void setTolerance(Scalar value) { params_.tolerance_ = value; } + /*! + * \brief Public wrapper for applying a Newton update to a solution vector. + * + * This keeps update_() protected while allowing external orchestrators + * (like flow-level nonlinear-system wrappers) to request a model-consistent + * update step. + */ + void applyUpdate(SolutionVector& nextSolution, + const SolutionVector& currentSolution, + const GlobalEqVector& solutionUpdate, + const GlobalEqVector& currentResidual) + { + asImp_().update_(nextSolution, currentSolution, solutionUpdate, currentResidual); + } + /*! * \brief Run the Newton method. * diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh index 89103ef70fd..6c6699afcca 100644 --- a/opm/models/ptflash/flashprimaryvariables.hh +++ b/opm/models/ptflash/flashprimaryvariables.hh @@ -90,6 +90,14 @@ public: using ParentType::operator=; //!< Import base class assignment operators. + template + void serializeOp(Serializer& serializer) + { + for (unsigned i = 0; i < this->size(); ++i) { + serializer((*this)[i]); + } + } + /*! * \copydoc ImmisciblePrimaryVariables::assignMassConservative */ diff --git a/opm/models/utils/parametersystem.cpp b/opm/models/utils/parametersystem.cpp index ff17aebb898..fa299fe2de1 100644 --- a/opm/models/utils/parametersystem.cpp +++ b/opm/models/utils/parametersystem.cpp @@ -424,6 +424,23 @@ void Hide_(const std::string& paramName) paramInfo.isHidden = true; } +bool HideIfRegistered_(const std::string& paramName) +{ + if (!MetaData::registrationOpen()) { + throw std::logic_error("Parameter '" + paramName + "' declared as hidden" + " when parameter registration was already closed."); + } + + auto paramInfoIt = MetaData::mutableRegistry().find(paramName); + if (paramInfoIt == MetaData::mutableRegistry().end()) { + return false; + } + + auto& paramInfo = paramInfoIt->second; + paramInfo.isHidden = true; + return true; +} + bool IsSet_(const std::string& paramName, bool errorIfNotRegistered) { if (errorIfNotRegistered) { diff --git a/opm/models/utils/parametersystem.hpp b/opm/models/utils/parametersystem.hpp index 1e4a4ddfc5a..ca71902efa9 100644 --- a/opm/models/utils/parametersystem.hpp +++ b/opm/models/utils/parametersystem.hpp @@ -78,6 +78,9 @@ ParamType Get_(const std::string& paramName, ParamType defaultValue, //! \brief Private implementation. void Hide_(const std::string& paramName); +//! \brief Private implementation. +bool HideIfRegistered_(const std::string& paramName); + //! \brief Private implementation. bool IsSet_(const std::string& paramName, bool errorIfNotRegistered); @@ -313,6 +316,17 @@ void Hide() detail::Hide_(detail::getParamName()); } +/*! + * \brief Indicate that a given parameter should not be mentioned in the help message if present. + * + * This is useful for parameters that are only registered for some models or output modules. + */ +template +bool HideIfRegistered() +{ + return detail::HideIfRegistered_(detail::getParamName()); +} + /*! * \brief Query whether parameter registration is open or not. * \return True if registration is open, false otherwise diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 41697958de0..b42f2d463eb 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -190,7 +190,7 @@ struct NupcolGroupRateTolerance { static constexpr Scalar value = 0.001; }; namespace Opm { -/// Solver parameters for the BlackoilModel. +/// Solver parameters for the NonlinearSystemBlackOilReservoir. template struct BlackoilModelParameters { diff --git a/opm/simulators/flow/FlowGenericProblem.cpp b/opm/simulators/flow/FlowGenericProblem.cpp index 45b64eab304..76e076c6010 100644 --- a/opm/simulators/flow/FlowGenericProblem.cpp +++ b/opm/simulators/flow/FlowGenericProblem.cpp @@ -29,6 +29,7 @@ #include #include +#include #include @@ -40,29 +41,81 @@ namespace Opm { +using CpLeafGridView = Dune::GridView>; + #define INSTANTIATE_TYPE(T) \ template class FlowGenericProblem< \ - Dune::GridView< \ - Dune::DefaultLeafGridViewTraits>, \ + CpLeafGridView, \ BlackOilFluidSystem>; +#define INSTANTIATE_COMP_TYPE(T, N, W) \ + template class FlowGenericProblem< \ + CpLeafGridView, \ + GenericOilGasWaterFluidSystem>; + +#define INSTANTIATE_COMP_VARIANTS(T, N) \ + INSTANTIATE_COMP_TYPE(T, N, false) \ + INSTANTIATE_COMP_TYPE(T, N, true) + +#define INSTANTIATE_FLOW_GENERIC_COMPONENTS(T) \ + INSTANTIATE_COMP_VARIANTS(T, 2) \ + INSTANTIATE_COMP_VARIANTS(T, 3) \ + INSTANTIATE_COMP_VARIANTS(T, 4) \ + INSTANTIATE_COMP_VARIANTS(T, 5) \ + INSTANTIATE_COMP_VARIANTS(T, 6) \ + INSTANTIATE_COMP_VARIANTS(T, 7) \ + INSTANTIATE_COMP_VARIANTS(T, 8) \ + INSTANTIATE_COMP_VARIANTS(T, 9) \ + INSTANTIATE_COMP_VARIANTS(T, 10) + INSTANTIATE_TYPE(double) +INSTANTIATE_FLOW_GENERIC_COMPONENTS(double) #if FLOW_INSTANTIATE_FLOAT INSTANTIATE_TYPE(float) +INSTANTIATE_FLOW_GENERIC_COMPONENTS(float) #endif #if HAVE_DUNE_FEM using GV = Dune::Fem::AdaptiveLeafGridPart; + +#define INSTANTIATE_DUNE_FEM_COMP_TYPE(T, N, W) \ +template class FlowGenericProblem>; + +#define INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, N) \ + INSTANTIATE_DUNE_FEM_COMP_TYPE(T, N, false) \ + INSTANTIATE_DUNE_FEM_COMP_TYPE(T, N, true) + +#define INSTANTIATE_DUNE_FEM_FLOW_GENERIC_COMPONENTS(T) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 2) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 3) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 4) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 5) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 6) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 7) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 8) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 9) \ + INSTANTIATE_DUNE_FEM_COMP_VARIANTS(T, 10) + template class FlowGenericProblem>; +INSTANTIATE_DUNE_FEM_FLOW_GENERIC_COMPONENTS(double) #if FLOW_INSTANTIATE_FLOAT template class FlowGenericProblem>; +INSTANTIATE_DUNE_FEM_FLOW_GENERIC_COMPONENTS(float) #endif +#undef INSTANTIATE_DUNE_FEM_FLOW_GENERIC_COMPONENTS +#undef INSTANTIATE_DUNE_FEM_COMP_VARIANTS +#undef INSTANTIATE_DUNE_FEM_COMP_TYPE + #endif // HAVE_DUNE_FEM +#undef INSTANTIATE_FLOW_GENERIC_COMPONENTS +#undef INSTANTIATE_COMP_VARIANTS + } // end namespace Opm diff --git a/opm/simulators/flow/FlowMain.hpp b/opm/simulators/flow/FlowMain.hpp index ff582e7d622..6939fd5c170 100644 --- a/opm/simulators/flow/FlowMain.hpp +++ b/opm/simulators/flow/FlowMain.hpp @@ -31,7 +31,7 @@ #include #include #include -#include +#include #include #if HAVE_DUNE_FEM @@ -61,7 +61,7 @@ namespace Opm { class Deck; - // The FlowMain class is the black-oil simulator. + // The FlowMain class is the standard fully implicit flow simulator. template class FlowMain { @@ -74,7 +74,7 @@ namespace Opm { using Scalar = GetPropType; using FluidSystem = GetPropType; - using Simulator = SimulatorFullyImplicitBlackoil; + using Simulator = SimulatorFullyImplicit; FlowMain(int argc, char **argv, bool output_cout, bool output_files ) : argc_{argc}, argv_{argv}, diff --git a/opm/simulators/flow/FlowProblemBlackoilProperties.hpp b/opm/simulators/flow/FlowProblemBlackoilProperties.hpp index 8f58a50fd6c..690efa41c92 100644 --- a/opm/simulators/flow/FlowProblemBlackoilProperties.hpp +++ b/opm/simulators/flow/FlowProblemBlackoilProperties.hpp @@ -37,7 +37,7 @@ #include #include #include -#include +#include #include @@ -60,7 +60,7 @@ struct FlowBaseProblemBlackoil { // Set the problem property template struct NonlinearSystem -{ using type = BlackoilModel; }; +{ using type = NonlinearSystemBlackOilReservoir; }; template diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index 85b9ef659ca..7087efe4e1e 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -424,6 +424,18 @@ class FlowProblemComp : public FlowProblem EclWriterType& eclWriter() { return *eclWriter_; } + void setSubStepReport(const SimulatorReportSingle& report) + { return eclWriter_->setSubStepReport(report); } + + void setSimulationReport(const SimulatorReport& report) + { return eclWriter_->setSimulationReport(report); } + + void finalizeOutput() + { + OPM_TIMEBLOCK(finalizeOutput); + eclWriter_.reset(); + } + // TODO: do we need this one? template void serializeOp(Serializer& serializer) diff --git a/opm/simulators/flow/FlowProblemCompProperties.hpp b/opm/simulators/flow/FlowProblemCompProperties.hpp index 17128bf019c..d671d368b6f 100644 --- a/opm/simulators/flow/FlowProblemCompProperties.hpp +++ b/opm/simulators/flow/FlowProblemCompProperties.hpp @@ -33,6 +33,7 @@ #include #include +#include #include #include @@ -51,6 +52,10 @@ struct FlowBaseProblemComp { }; } +template +struct NonlinearSystem +{ using type = NonlinearSystemCompositional; }; + // Set the problem property template struct Problem diff --git a/opm/simulators/flow/FlowUtils.cpp b/opm/simulators/flow/FlowUtils.cpp index 53de3495bdb..86d3e0b43a9 100644 --- a/opm/simulators/flow/FlowUtils.cpp +++ b/opm/simulators/flow/FlowUtils.cpp @@ -217,49 +217,49 @@ void hideUnusedParameters() Parameters::Hide(); // hide all vtk related it is not currently possible to do this dependet on if the vtk writing is used //if(not(Parameters::Get())){ - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); //} - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); - - Parameters::Hide(); - Parameters::Hide(); - Parameters::Hide(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); + Parameters::HideIfRegistered(); // hide average density option - Parameters::Hide(); + Parameters::HideIfRegistered(); } namespace { diff --git a/opm/simulators/flow/NonlinearSolver.cpp b/opm/simulators/flow/NonlinearSolver.cpp index 8203f01ac51..c2a194c30ea 100644 --- a/opm/simulators/flow/NonlinearSolver.cpp +++ b/opm/simulators/flow/NonlinearSolver.cpp @@ -120,7 +120,12 @@ using BV = Dune::BlockVector>; INSTANTIATE(T,3) \ INSTANTIATE(T,4) \ INSTANTIATE(T,5) \ - INSTANTIATE(T,6) + INSTANTIATE(T,6) \ + INSTANTIATE(T,7) \ + INSTANTIATE(T,8) \ + INSTANTIATE(T,9) \ + INSTANTIATE(T,10) \ + INSTANTIATE(T,11) INSTANTIATE_TYPE(double) diff --git a/opm/simulators/flow/NonlinearSystem.hpp b/opm/simulators/flow/NonlinearSystem.hpp new file mode 100644 index 00000000000..743b298c77a --- /dev/null +++ b/opm/simulators/flow/NonlinearSystem.hpp @@ -0,0 +1,181 @@ +#ifndef OPM_FLOW_NONLINEAR_SYSTEM_HEADER_INCLUDED +#define OPM_FLOW_NONLINEAR_SYSTEM_HEADER_INCLUDED +/* + Copyright 2026, 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 + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +template +class NonlinearSystem +{ +public: + // --------- Types --------- + using Simulator = GetPropType; + using Grid = GetPropType; + using FluidSystem = GetPropType; + using Indices = GetPropType; + using Scalar = GetPropType; + using GlobalEqVector = GetPropType; + using ModelParameters = BlackoilModelParameters; + using WellModel = GetPropType; + using ComponentName = ::Opm::ComponentName; + + // --------- Public methods --------- + virtual ~NonlinearSystem() = default; + + bool isParallel() const + { return grid_.comm().size() > 1; } + + const Simulator& simulator() const + { return simulator_; } + + Simulator& simulator() + { return simulator_; } + + bool terminalOutputEnabled() const + { return terminal_output_; } + + int numPhases() const + { return Indices::numPhases; } + + const SimulatorReportSingle& failureReport() const + { return failureReport_; } + + const std::vector& stepReports() const + { return convergence_reports_; } + + const ComponentName& compNames() const + { return compNames_; } + + const ModelParameters& param() const + { return param_; } + + WellModel& wellModel() + { return well_model_; } + + const WellModel& wellModel() const + { return well_model_; } + + void beginReportStep() + { simulator_.problem().beginEpisode(); } + + void endReportStep() + { simulator_.problem().endEpisode(); } + + SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface& timer); + + void updateTUNING(const Tuning& tuning); + + void updateTUNINGDP(const TuningDp& tuning_dp); + + void updateSolution(const GlobalEqVector& dx); + + template + void addReservoirConvergenceMetrics(ConvergenceReport& report, + const int componentIdx, + const std::string_view componentName, + const std::span residuals, + const std::span types, + const std::span tolerances, + const Scalar maxResidualAllowed, + LogFailure&& logFailure) const; + +protected: + explicit NonlinearSystem(Simulator& simulator, + const ModelParameters& param, + WellModel& wellModel, + const bool terminal_output); + + virtual void initialLinearization(SimulatorReportSingle& report, + int minIter, + int maxIter, + const SimulatorTimerInterface& timer); + + virtual bool shouldStoreSolutionUpdate() const + { return false; } + + virtual void prepareSolutionUpdate() + {} + + virtual void storeSolutionUpdate(const GlobalEqVector&) + {} + + SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer); + + template + SimulatorReportSingle assembleReservoir(WellModelType& wellModel); + + template + void applyTUNING(ModelParametersType& param, + const Tuning& tuning); + + template + void applyTUNINGDP(ModelParametersType& param, + const TuningDp& tuning_dp); + + template + std::tuple + convergenceReduction(Parallel::Communication comm, + const ValueType primaryVolumeLocal, + const ValueType secondaryVolumeLocal, + std::vector& sumValues, + std::vector& maxValues, + std::vector& averagedValues); + + void popLastStepReport() + { convergence_reports_.back().report.pop_back(); } + + // --------- Data members --------- + Simulator& simulator_; + const Grid& grid_; + bool terminal_output_; + ModelParameters param_; + WellModel& well_model_; + SimulatorReportSingle failureReport_; + std::vector convergence_reports_; + ComponentName compNames_{}; + std::vector> residual_norms_history_; + Scalar current_relaxation_; + GlobalEqVector dx_old_; +}; + +} // namespace Opm + +#include + +#endif // OPM_FLOW_NONLINEAR_SYSTEM_HEADER_INCLUDED diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoir.hpp similarity index 83% rename from opm/simulators/flow/BlackoilModel.hpp rename to opm/simulators/flow/NonlinearSystemBlackOilReservoir.hpp index e13aed816d0..9dca9de2218 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoir.hpp @@ -21,13 +21,15 @@ along with OPM. If not, see . */ -#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED -#define OPM_BLACKOILMODEL_HEADER_INCLUDED +#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_HEADER_INCLUDED +#define OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_HEADER_INCLUDED #include +#include + #include -#include +#include #include #include #include @@ -57,12 +59,13 @@ namespace Opm { /// uses an industry-standard TPFA discretization with per-phase /// upwind weighting of mobilities. template -class BlackoilModel +class NonlinearSystemBlackOilReservoir : public NonlinearSystem { public: // --------- Types and enums --------- - using Simulator = GetPropType; - using Grid = GetPropType; + using ParentType = NonlinearSystem; + using Simulator = typename ParentType::Simulator; + using Grid = typename ParentType::Grid; using ElementContext = GetPropType; using IntensiveQuantities = GetPropType; using SparseMatrixAdapter = GetPropType; @@ -73,6 +76,7 @@ class BlackoilModel using MaterialLaw = GetPropType; using MaterialLawParams = GetPropType; using Scalar = GetPropType; + using GlobalEqVector = typename ParentType::GlobalEqVector; using ModelParameters = BlackoilModelParameters; static constexpr bool enableSaltPrecipitation = getPropValue(); @@ -139,14 +143,11 @@ class BlackoilModel /// \param[in] param parameters /// \param[in] well_model Reference to well model /// \param[in] terminal_output request output to cout/cerr - BlackoilModel(Simulator& simulator, + NonlinearSystemBlackOilReservoir(Simulator& simulator, const ModelParameters& param, BlackoilWellModel& well_model, const bool terminal_output); - bool isParallel() const - { return grid_.comm().size() > 1; } - const EclipseState& eclState() const { return simulator_.vanguard().eclState(); } @@ -157,7 +158,7 @@ class BlackoilModel void initialLinearization(SimulatorReportSingle& report, const int minIter, const int maxIter, - const SimulatorTimerInterface& timer); + const SimulatorTimerInterface& timer) override; /// Called once per nonlinear iteration. /// This model will perform a Newton-Raphson update, changing reservoir_state @@ -174,9 +175,6 @@ class BlackoilModel SimulatorReportSingle nonlinearIterationNewton(const SimulatorTimerInterface& timer, NonlinearSolverType& nonlinear_solver); - /// Assemble the residual and Jacobian of the nonlinear system. - SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface& timer); - // compute the "relative" change of the solution between time steps Scalar relativeChange() const; @@ -192,18 +190,12 @@ class BlackoilModel /// r is the residual. void solveJacobianSystem(BVector& x); - /// Apply an update to the primary variables. - void updateSolution(const BVector& dx); - /// Get solution update vector as a PrimaryVariable - void prepareStoringSolutionUpdate(); - void storeSolutionUpdate(const BVector& dx); + bool shouldStoreSolutionUpdate() const override; + void prepareSolutionUpdate() override; + void storeSolutionUpdate(const GlobalEqVector& dx) override; MaxSolutionUpdateData getMaxSolutionUpdate(const std::vector& ixCells); - /// Return true if output to cout is wanted. - bool terminalOutputEnabled() const - { return terminal_output_; } - std::tuple convergenceReduction(Parallel::Communication comm, const Scalar pvSumLocal, @@ -234,9 +226,6 @@ class BlackoilModel const double tol_cnv, const double tol_cnv_energy); - void updateTUNING(const Tuning& tuning); - void updateTUNINGDP(const TuningDp& tuning_dp); - ConvergenceReport getReservoirConvergence(const double reportTime, const double dt, @@ -254,10 +243,6 @@ class BlackoilModel const int maxIter, std::vector& residual_norms); - /// The number of active fluid phases in the model. - int numPhases() const - { return Indices::numPhases; } - /// Wrapper required due to not following generic API template std::vector > @@ -268,16 +253,6 @@ class BlackoilModel std::vector > computeFluidInPlace(const std::vector& /*fipnum*/) const; - const Simulator& simulator() const - { return simulator_; } - - Simulator& simulator() - { return simulator_; } - - /// return the statistics if the nonlinearIteration() method failed - const SimulatorReportSingle& failureReport() const - { return failureReport_; } - /// return the statistics of local solves accumulated for this rank const SimulatorReport& localAccumulatedReports() const; @@ -287,33 +262,15 @@ class BlackoilModel /// Write the number of nonlinear iterations per cell to a file in ResInsight compatible format void writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const; - const std::vector& stepReports() const - { return convergence_reports_; } - /// Remove the last convergence report entry and residual norms history entry. void popLastConvergenceReport() { - convergence_reports_.back().report.pop_back(); - residual_norms_history_.pop_back(); + this->popLastStepReport(); + this->residual_norms_history_.pop_back(); } void writePartitions(const std::filesystem::path& odir) const; - /// return the StandardWells object - BlackoilWellModel& - wellModel() - { return well_model_; } - - const BlackoilWellModel& - wellModel() const - { return well_model_; } - - void beginReportStep() - { simulator_.problem().beginEpisode(); } - - void endReportStep() - { simulator_.problem().endEpisode(); } - template void getMaxCoeff(const unsigned cell_idx, const IntensiveQuantities& intQuants, @@ -325,22 +282,25 @@ class BlackoilModel std::vector& maxCoeff, std::vector& maxCoeffCell); - //! \brief Returns const reference to model parameters. - const ModelParameters& param() const - { return param_; } - - //! \brief Returns const reference to component names. - const ComponentName& compNames() const - { return compNames_; } - //! \brief Returns true if an NLDD solver exists bool hasNlddSolver() const { return nlddSolver_ != nullptr; } protected: + using ParentType::compNames_; + using ParentType::convergence_reports_; + using ParentType::current_relaxation_; + using ParentType::dx_old_; + using ParentType::failureReport_; + using ParentType::grid_; + using ParentType::param_; + using ParentType::popLastStepReport; + using ParentType::residual_norms_history_; + using ParentType::simulator_; + using ParentType::terminal_output_; + using ParentType::well_model_; + // --------- Data members --------- - Simulator& simulator_; - const Grid& grid_; static constexpr bool has_solvent_ = getPropValue(); static constexpr bool has_extbo_ = getPropValue(); static constexpr bool has_polymer_ = getPropValue(); @@ -351,27 +311,12 @@ class BlackoilModel static constexpr bool has_bioeffects_ = getPropValue(); static constexpr bool has_micp_ = Indices::enableMICP; - ModelParameters param_; - SimulatorReportSingle failureReport_; - - // Well Model - BlackoilWellModel& well_model_; - - /// \brief Whether we print something to std::cout - bool terminal_output_; /// \brief The number of cells of the global grid. long int global_nc_; - std::vector> residual_norms_history_; - Scalar current_relaxation_; - BVector dx_old_; - SolutionVector solUpd_; - std::vector convergence_reports_; - ComponentName compNames_{}; - - std::unique_ptr> nlddSolver_; //!< Non-linear DD solver + std::unique_ptr> nlddSolver_; //!< Non-linear DD solver BlackoilModelConvergenceMonitor conv_monitor_; private: @@ -385,6 +330,6 @@ class BlackoilModel } // namespace Opm -#include +#include -#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED +#endif // OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_HEADER_INCLUDED diff --git a/opm/simulators/flow/BlackoilModelTPSA.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp similarity index 92% rename from opm/simulators/flow/BlackoilModelTPSA.hpp rename to opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp index 620761a43ab..c4d33969828 100644 --- a/opm/simulators/flow/BlackoilModelTPSA.hpp +++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp @@ -22,13 +22,13 @@ module for the precise wording of the license and the list of copyright holders. */ -#ifndef BLACK_OIL_MODEL_TPSA_HPP -#define BLACK_OIL_MODEL_TPSA_HPP +#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_TPSA_HPP +#define OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_TPSA_HPP #include #include -#include +#include #include #include @@ -39,12 +39,12 @@ namespace Opm { /*! -* \brief Black oil model for coupling Flow simulations with TPSA geomechanics +* \brief Nonlinear system for coupling Flow simulations with TPSA geomechanics */ template -class BlackoilModelTPSA : public BlackoilModel +class NonlinearSystemBlackOilReservoirTPSA : public NonlinearSystemBlackOilReservoir { - using ParentType = BlackoilModel; + using ParentType = NonlinearSystemBlackOilReservoir; using Scalar = GetPropType; using Simulator = GetPropType; @@ -60,7 +60,7 @@ class BlackoilModelTPSA : public BlackoilModel * \param well_model Refenerence to well model * \param terminal_output Bool for terminal output */ - explicit BlackoilModelTPSA(Simulator& simulator, + explicit NonlinearSystemBlackOilReservoirTPSA(Simulator& simulator, const ModelParameters& param, BlackoilWellModel& well_model, const bool terminal_output) @@ -223,8 +223,8 @@ class BlackoilModelTPSA : public BlackoilModel private: int seqIter_{0}; -}; // class BlackoilModelTPSA +}; // class NonlinearSystemBlackOilReservoirTPSA } // namespace Opm -#endif +#endif // OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_TPSA_HPP diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp similarity index 78% rename from opm/simulators/flow/BlackoilModel_impl.hpp rename to opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp index b680e76ea11..6f59afb3549 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp @@ -21,12 +21,12 @@ along with OPM. If not, see . */ -#ifndef OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED -#define OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED +#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED +#define OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED -#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED +#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_HEADER_INCLUDED #include -#include +#include #endif #include @@ -57,9 +57,9 @@ namespace { template std::string_view - make_string(const typename Opm::BlackoilModel::DebugFlags f) + make_string(const typename Opm::NonlinearSystemBlackOilReservoir::DebugFlags f) { - using F = typename Opm::BlackoilModel::DebugFlags; + using F = typename Opm::NonlinearSystemBlackOilReservoir::DebugFlags; switch (f) { case F::STRICT: return "Strict"; @@ -74,19 +74,13 @@ namespace { namespace Opm { template -BlackoilModel:: -BlackoilModel(Simulator& simulator, +NonlinearSystemBlackOilReservoir:: +NonlinearSystemBlackOilReservoir(Simulator& simulator, const ModelParameters& param, BlackoilWellModel& well_model, const bool terminal_output) - : simulator_(simulator) - , grid_(simulator_.vanguard().grid()) - , param_( param ) - , well_model_ (well_model) - , terminal_output_ (terminal_output) - , current_relaxation_(1.0) - , dx_old_(simulator_.model().numGridDof()) - , conv_monitor_(param_.monitor_params_) + : ParentType(simulator, param, well_model, terminal_output) + , conv_monitor_(param.monitor_params_) { // compute global sum of number of cells global_nc_ = detail::countGlobalCells(grid_); @@ -96,7 +90,7 @@ BlackoilModel(Simulator& simulator, if (terminal_output) { OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd)."); } - nlddSolver_ = std::make_unique>(*this); + nlddSolver_ = std::make_unique>(*this); } else if (param_.nonlinear_solver_ == "newton") { if (terminal_output) { OpmLog::info("Using Newton nonlinear solver."); @@ -109,36 +103,14 @@ BlackoilModel(Simulator& simulator, template SimulatorReportSingle -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: prepareStep(const SimulatorTimerInterface& timer) { OPM_TIMEFUNCTION(); - SimulatorReportSingle report; + auto report = ParentType::prepareStep(timer); + Dune::Timer perfTimer; perfTimer.start(); - // update the solution variables in the model - int lastStepFailed = timer.lastStepFailed(); - if (grid_.comm().size() > 1 && grid_.comm().max(lastStepFailed) != grid_.comm().min(lastStepFailed)) { - OPM_THROW(std::runtime_error, "Misalignment of the parallel simulation run in prepareStep " + - "- the previous step succeeded on some ranks but failed on others."); - } - if (lastStepFailed) { - simulator_.model().updateFailed(); - } - else { - simulator_.model().advanceTimeLevel(); - } - - // Set the timestep size and episode index for the model explicitly. - // The model needs to know the report step/episode index because of - // timing dependent data despite the fact that flow uses its own time stepper. - // (The length of the episode does not matter, though.) - simulator_.setTime(timer.simulationTimeElapsed()); - simulator_.setTimeStepSize(timer.currentStepLength()); - - simulator_.problem().resetIterationForNewTimestep(); - - simulator_.problem().beginTimeStep(); unsigned numDof = simulator_.model().numGridDof(); wasSwitched_.resize(numDof); @@ -180,35 +152,20 @@ prepareStep(const SimulatorTimerInterface& timer) template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: initialLinearization(SimulatorReportSingle& report, const int minIter, const int maxIter, const SimulatorTimerInterface& timer) { - // ----------- Set up reports and timer ----------- - failureReport_ = SimulatorReportSingle(); - Dune::Timer perfTimer; - - perfTimer.start(); - report.total_linearizations = 1; - - // ----------- Assemble ----------- - try { - report += assembleReservoir(timer); - report.assemble_time += perfTimer.stop(); - // Mark timestep initialized after assembleReservoir(), because the well model's - // assemble() also uses needsTimestepInit() to trigger prepareTimeStep(). - simulator_.problem().markTimestepInitialized(); - } - catch (...) { - report.assemble_time += perfTimer.stop(); - failureReport_ += report; - throw; // continue throwing the stick - } + ParentType::initialLinearization(report, + minIter, + maxIter, + timer); // ----------- Check if converged ----------- std::vector residual_norms; + Dune::Timer perfTimer; perfTimer.reset(); perfTimer.start(); // the step is not considered converged until at least minIter iterations is done @@ -242,48 +199,46 @@ initialLinearization(SimulatorReportSingle& report, template template SimulatorReportSingle -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: nonlinearIteration(const SimulatorTimerInterface& timer, NonlinearSolverType& nonlinear_solver) { // Model-level timestep initialization (once per timestep). // markTimestepInitialized() is called later in initialLinearization(), // after assembleReservoir() has triggered the well model's prepareTimeStep(). - if (simulator_.problem().iterationContext().needsTimestepInit()) { + if (this->simulator_.problem().iterationContext().needsTimestepInit()) { residual_norms_history_.clear(); conv_monitor_.reset(); current_relaxation_ = 1.0; dx_old_ = 0.0; - convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}}); - convergence_reports_.back().report.reserve(11); + this->convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}}); + this->convergence_reports_.back().report.reserve(11); } SimulatorReportSingle result; - if (this->param_.nonlinear_solver_ != "nldd") - { + if (this->param_.nonlinear_solver_ != "nldd") { result = this->nonlinearIterationNewton(timer, nonlinear_solver); } else { result = this->nlddSolver_->nonlinearIterationNldd(timer, nonlinear_solver); } - auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv(); - rst_conv.update(simulator_.model().linearizer().residual()); + auto& rst_conv = this->simulator_.problem().eclWriter().mutableOutputModule().getConv(); + rst_conv.update(this->simulator_.model().linearizer().residual()); - simulator_.problem().advanceIteration(); + this->simulator_.problem().advanceIteration(); return result; } template template SimulatorReportSingle -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: nonlinearIterationNewton(const SimulatorTimerInterface& timer, NonlinearSolverType& nonlinear_solver) { OPM_TIMEFUNCTION(); - // ----------- Set up reports and timer ----------- SimulatorReportSingle report; Dune::Timer perfTimer; @@ -292,26 +247,19 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, this->param_.newton_max_iter_, timer); - // ----------- If not converged, solve linear system and do Newton update ----------- if (!report.converged) { perfTimer.reset(); perfTimer.start(); report.total_newton_iterations = 1; - // Compute the nonlinear update. - unsigned nc = simulator_.model().numGridDof(); + const unsigned nc = this->simulator_.model().numGridDof(); BVector x(nc); - // Solve the linear system. linear_solve_setup_time_ = 0.0; try { - // Apply the Schur complement of the well model to - // the reservoir linearized equations. - // Note that linearize may throw for MSwells. - wellModel().linearize(simulator().model().linearizer().jacobian(), - simulator().model().linearizer().residual()); + this->wellModel().linearize(this->simulator().model().linearizer().jacobian(), + this->simulator().model().linearizer().residual()); - // ---- Solve linear system ---- solveJacobianSystem(x); report.linear_solve_setup_time += linear_solve_setup_time_; @@ -323,20 +271,16 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, report.linear_solve_time += perfTimer.stop(); report.total_linear_iterations += linearIterationsLastSolve(); - failureReport_ += report; - throw; // re-throw up + this->failureReport_ += report; + throw; } perfTimer.reset(); perfTimer.start(); - // handling well state update before oscillation treatment is a decision based - // on observation to avoid some big performance degeneration under some circumstances. - // there is no theorectical explanation which way is better for sure. - wellModel().postSolve(x); + this->wellModel().postSolve(x); if (param_.use_update_stabilization_) { - // Stabilize the nonlinear update. bool isOscillate = false; bool isStagnate = false; nonlinear_solver.detectOscillations(residual_norms_history_, @@ -345,22 +289,16 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, isStagnate); if (isOscillate) { current_relaxation_ -= nonlinear_solver.relaxIncrement(); - current_relaxation_ = std::max(current_relaxation_, - nonlinear_solver.relaxMax()); - if (terminalOutputEnabled()) { - std::string msg = " Oscillating behavior detected: Relaxation set to " - + std::to_string(current_relaxation_); - OpmLog::info(msg); + current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); + if (this->terminalOutputEnabled()) { + OpmLog::info(" Oscillating behavior detected: Relaxation set to " + + std::to_string(current_relaxation_)); } } nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_); } - // ---- Newton update ---- - // Apply the update, with considering model-dependent limitations and - // chopping of the update. - updateSolution(x); - + this->updateSolution(x); report.update_time += perfTimer.stop(); } @@ -368,20 +306,8 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, } template -SimulatorReportSingle -BlackoilModel:: -assembleReservoir(const SimulatorTimerInterface& /* timer */) -{ - // -------- Mass balance equations -------- - simulator_.problem().beginIteration(); - simulator_.model().linearizer().linearizeDomain(); - simulator_.problem().endIteration(); - return wellModel().lastReport(); -} - -template -typename BlackoilModel::Scalar -BlackoilModel:: +typename NonlinearSystemBlackOilReservoir::Scalar +NonlinearSystemBlackOilReservoir:: relativeChange() const { Scalar resultDelta = 0.0; @@ -468,7 +394,7 @@ relativeChange() const template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: solveJacobianSystem(BVector& x) { auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix(); @@ -488,7 +414,6 @@ solveJacobianSystem(BVector& x) x = 0.0; std::vector x_trial(numSolvers, x); for (int solver = 0; solver < numSolvers; ++solver) { - BVector x0(x); linSolver.setActiveSolver(solver); perfTimer.start(); linSolver.prepare(jacobian, residual); @@ -512,7 +437,6 @@ solveJacobianSystem(BVector& x) linSolver.setActiveSolver(fastest_solver); } else { - // set initial guess x = 0.0; Dune::Timer perfTimer; @@ -529,44 +453,18 @@ solveJacobianSystem(BVector& x) } template -void -BlackoilModel:: -updateSolution(const BVector& dx) +bool +NonlinearSystemBlackOilReservoir:: +shouldStoreSolutionUpdate() const { - OPM_TIMEBLOCK(updateSolution); - // Prepare to store solution update for convergence check - if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 - || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) { - prepareStoringSolutionUpdate(); - } - - auto& newtonMethod = simulator_.model().newtonMethod(); - SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0); - - newtonMethod.update_(/*nextSolution=*/solution, - /*curSolution=*/solution, - /*update=*/dx, - /*resid=*/dx); // the update routines of the black - // oil model do not care about the - // residual - - // if the solution is updated, the intensive quantities need to be recalculated - { - OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities); - simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); - } - - // Store solution update - if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 - || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) { - storeSolutionUpdate(dx); - } + return this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0; } template void -BlackoilModel:: -prepareStoringSolutionUpdate() +NonlinearSystemBlackOilReservoir:: +prepareSolutionUpdate() { // Init. solution update vector unsigned nc = simulator_.model().numGridDof(); @@ -586,8 +484,8 @@ prepareStoringSolutionUpdate() template void -BlackoilModel:: -storeSolutionUpdate(const BVector& dx) +NonlinearSystemBlackOilReservoir:: +storeSolutionUpdate(const GlobalEqVector& dx) { const auto& elemMapper = simulator_.model().elementMapper(); const auto& gridView = simulator_.gridView(); @@ -604,8 +502,8 @@ storeSolutionUpdate(const BVector& dx) } template -typename BlackoilModel::MaxSolutionUpdateData -BlackoilModel:: +typename NonlinearSystemBlackOilReservoir::MaxSolutionUpdateData +NonlinearSystemBlackOilReservoir:: getMaxSolutionUpdate(const std::vector& ixCells) { static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0; @@ -655,9 +553,9 @@ getMaxSolutionUpdate(const std::vector& ixCells) } template -std::tuple::Scalar, - typename BlackoilModel::Scalar> -BlackoilModel:: +std::tuple::Scalar, + typename NonlinearSystemBlackOilReservoir::Scalar> +NonlinearSystemBlackOilReservoir:: convergenceReduction(Parallel::Communication comm, const Scalar pvSumLocal, const Scalar numAquiferPvSumLocal, @@ -665,59 +563,18 @@ convergenceReduction(Parallel::Communication comm, std::vector< Scalar >& maxCoeff, std::vector< Scalar >& B_avg) { - OPM_TIMEBLOCK(convergenceReduction); - // Compute total pore volume (use only owned entries) - Scalar pvSum = pvSumLocal; - Scalar numAquiferPvSum = numAquiferPvSumLocal; - - if (comm.size() > 1) { - // global reduction - std::vector< Scalar > sumBuffer; - std::vector< Scalar > maxBuffer; - const int numComp = B_avg.size(); - sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum - maxBuffer.reserve( numComp ); - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - sumBuffer.push_back(B_avg[compIdx]); - sumBuffer.push_back(R_sum[compIdx]); - maxBuffer.push_back(maxCoeff[ compIdx]); - } - - // Compute total pore volume - sumBuffer.push_back( pvSum ); - sumBuffer.push_back( numAquiferPvSum ); - - // compute global sum - comm.sum( sumBuffer.data(), sumBuffer.size() ); - - // compute global max - comm.max( maxBuffer.data(), maxBuffer.size() ); - - // restore values to local variables - for (int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx) { - B_avg[compIdx] = sumBuffer[buffIdx]; - ++buffIdx; - - R_sum[compIdx] = sumBuffer[buffIdx]; - } - - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - maxCoeff[compIdx] = maxBuffer[compIdx]; - } - - // restore global pore volume - pvSum = sumBuffer[sumBuffer.size()-2]; - numAquiferPvSum = sumBuffer.back(); - } - - // return global pore volume - return {pvSum, numAquiferPvSum}; + return ParentType::convergenceReduction(comm, + pvSumLocal, + numAquiferPvSumLocal, + R_sum, + maxCoeff, + B_avg); } template -std::pair::Scalar, - typename BlackoilModel::Scalar> -BlackoilModel:: +std::pair::Scalar, + typename NonlinearSystemBlackOilReservoir::Scalar> +NonlinearSystemBlackOilReservoir:: localConvergenceData(std::vector& R_sum, std::vector& maxCoeff, std::vector& B_avg, @@ -732,7 +589,7 @@ localConvergenceData(std::vector& R_sum, const auto& residual = simulator_.model().linearizer().residual(); ElementContext elemCtx(simulator_); - const auto& gridView = simulator().gridView(); + const auto& gridView = this->simulator().gridView(); IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { @@ -755,7 +612,7 @@ localConvergenceData(std::vector& R_sum, B_avg, R_sum, maxCoeff, maxCoeffCell); } - OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::localConvergenceData() failed: ", grid_.comm()); + OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::localConvergenceData() failed: ", grid_.comm()); // compute local average in terms of global number of elements const int bSize = B_avg.size(); @@ -767,8 +624,8 @@ localConvergenceData(std::vector& R_sum, } template -typename BlackoilModel::CnvPvSplitData -BlackoilModel:: +typename NonlinearSystemBlackOilReservoir::CnvPvSplitData +NonlinearSystemBlackOilReservoir:: characteriseCnvPvSplit(const std::vector& B_avg, const double dt) { OPM_TIMEBLOCK(computeCnvErrorPv); @@ -838,7 +695,7 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) } } - OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ", + OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::characteriseCnvPvSplit() failed: ", this->grid_.comm()); this->grid_.comm().sum(splitPV .data(), splitPV .size()); @@ -847,34 +704,9 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) return { cnvPvSplit, ixCells }; } -template -void -BlackoilModel:: -updateTUNING(const Tuning& tuning) -{ - this->param_.tolerance_cnv_ = tuning.TRGCNV; - this->param_.tolerance_cnv_relaxed_ = tuning.XXXCNV; - this->param_.tolerance_mb_ = tuning.TRGMBE; - this->param_.tolerance_mb_relaxed_ = tuning.XXXMBE; - this->param_.newton_max_iter_ = tuning.NEWTMX; - this->param_.newton_min_iter_ = tuning.NEWTMN; -} - -template -void -BlackoilModel:: -updateTUNINGDP(const TuningDp& tuning_dp) -{ - // NOTE: If TUNINGDP item is _not_ set it should be 0.0 - this->param_.tolerance_max_dp_ = tuning_dp.TRGDDP; - this->param_.tolerance_max_ds_ = tuning_dp.TRGDDS; - this->param_.tolerance_max_drs_ = tuning_dp.TRGDDRS; - this->param_.tolerance_max_drv_ = tuning_dp.TRGDDRV; -} - template ConvergenceReport -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: getReservoirConvergence(const double reportTime, const double dt, const int maxIter, @@ -1012,31 +844,20 @@ getReservoirConvergence(const double reportTime, tol[1] = tol_cnv_energy; } - for (int ii : {0, 1}) { - if (std::isnan(res[ii])) { - report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); - if (this->terminal_output_) { - OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation."); - } - } - else if (res[ii] > maxResidualAllowed()) { - report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); - if (this->terminal_output_) { - OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation."); - } - } - else if (res[ii] < 0.0) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + this->addReservoirConvergenceMetrics( + report, + compIdx, + this->compNames_.name(compIdx), + std::span{res}, + std::span{types}, + std::span{tol}, + maxResidualAllowed(), + [this](const std::string& message) + { if (this->terminal_output_) { - OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation."); + OpmLog::debug(message); } - } - else if (res[ii] > tol[ii]) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); - } - - report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]); - } + }); } // Compute the Newton convergence per cell. @@ -1128,7 +949,7 @@ getReservoirConvergence(const double reportTime, template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: convergencePerCell(const std::vector& B_avg, const double dt, const double tol_cnv, @@ -1167,14 +988,14 @@ convergencePerCell(const std::vector& B_avg, } ++idx; } - OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::convergencePerCell() failed: ", + OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::convergencePerCell() failed: ", this->grid_.comm()); rst_conv.updateNewton(convNewt); } template ConvergenceReport -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: getConvergence(const SimulatorTimerInterface& timer, const int maxIter, std::vector& residual_norms) @@ -1187,8 +1008,8 @@ getConvergence(const SimulatorTimerInterface& timer, maxIter, B_avg, residual_norms); { OPM_TIMEBLOCK(getWellConvergence); - report += wellModel().getWellConvergence(B_avg, - /*checkWellGroupControlsAndNetwork*/report.converged()); + report += this->wellModel().getWellConvergence(B_avg, + /*checkWellGroupControlsAndNetwork*/report.converged()); } conv_monitor_.checkPenaltyCard(report, simulator_.problem().iterationContext().iteration()); @@ -1197,8 +1018,8 @@ getConvergence(const SimulatorTimerInterface& timer, } template -std::vector::Scalar> > -BlackoilModel:: +std::vector::Scalar> > +NonlinearSystemBlackOilReservoir:: computeFluidInPlace(const std::vector& /*fipnum*/) const { OPM_TIMEBLOCK(computeFluidInPlace); @@ -1210,7 +1031,7 @@ computeFluidInPlace(const std::vector& /*fipnum*/) const template const SimulatorReport& -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: localAccumulatedReports() const { if (!hasNlddSolver()) { @@ -1221,7 +1042,7 @@ localAccumulatedReports() const template const std::vector& -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: domainAccumulatedReports() const { if (!nlddSolver_) @@ -1231,7 +1052,7 @@ domainAccumulatedReports() const template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const { if (hasNlddSolver()) { @@ -1241,7 +1062,7 @@ writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: writePartitions(const std::filesystem::path& odir) const { if (hasNlddSolver()) { @@ -1268,7 +1089,7 @@ writePartitions(const std::filesystem::path& odir) const template template void -BlackoilModel:: +NonlinearSystemBlackOilReservoir:: getMaxCoeff(const unsigned cell_idx, const IntensiveQuantities& intQuants, const FluidState& fs, @@ -1390,4 +1211,4 @@ getMaxCoeff(const unsigned cell_idx, } // namespace Opm -#endif // OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED +#endif // OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/flow/NonlinearSystemCompositional.hpp b/opm/simulators/flow/NonlinearSystemCompositional.hpp new file mode 100644 index 00000000000..2c883278d22 --- /dev/null +++ b/opm/simulators/flow/NonlinearSystemCompositional.hpp @@ -0,0 +1,122 @@ +/* + Copyright 2026, 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 . +*/ + +#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED +#define OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +namespace Opm { + +template +class NonlinearSystemCompositional : public NonlinearSystem +{ +public: + using ParentType = NonlinearSystem; + using Simulator = typename ParentType::Simulator; + using Grid = typename ParentType::Grid; + using FluidSystem = typename ParentType::FluidSystem; + using Indices = typename ParentType::Indices; + using Scalar = typename ParentType::Scalar; + using ComponentName = typename ParentType::ComponentName; + using SparseMatrixAdapter = GetPropType; + using ModelParameters = BlackoilModelParameters; + + static constexpr int numEq = Indices::numEq; + + using VectorBlockType = Dune::FieldVector; + using BVector = Dune::BlockVector; + + NonlinearSystemCompositional(Simulator& simulator, + const ModelParameters& param, + CompWellModel& wellModel, + bool terminalOutput); + + SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer); + + void initialLinearization(SimulatorReportSingle& report, + int minIter, + int maxIter, + const SimulatorTimerInterface& timer) override; + + template + SimulatorReportSingle nonlinearIteration(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinearSolver); + + template + SimulatorReportSingle nonlinearIterationNewton(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinearSolver); + + Scalar relativeChange() const; + + int linearIterationsLastSolve() const + { return this->simulator_.model().newtonMethod().linearSolver().iterations(); } + + void solveJacobianSystem(BVector& x); + + bool hasNlddSolver() const + { return false; } + + const SimulatorReport& localAccumulatedReports() const + { + static const SimulatorReport emptyReport{}; + return emptyReport; + } + + const std::vector& domainAccumulatedReports() const + { + static const std::vector emptyReports{}; + return emptyReports; + } + + void writeNonlinearIterationsPerCell(const std::filesystem::path&) const {} + + template + std::vector> computeFluidInPlace(const T&, const std::vector& fipnum) const + { return computeFluidInPlace(fipnum); } + + std::vector> computeFluidInPlace(const std::vector&) const + { return {}; } + + void writePartitions(const std::filesystem::path&) const {} + +private: + std::vector reservoirResidualMetrics() const; + + double linear_solve_setup_time_ = 0.0; +}; + +} // namespace Opm + +#include + +#endif // OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED diff --git a/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp b/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp new file mode 100644 index 00000000000..db941b124db --- /dev/null +++ b/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp @@ -0,0 +1,335 @@ +/* + Copyright 2026, 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 . +*/ + +#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_IMPL_HEADER_INCLUDED +#define OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_IMPL_HEADER_INCLUDED + +#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED +#include +#include +#endif + +#include + +#include +#include + +#include +#include +#include + +namespace Opm { + +template +NonlinearSystemCompositional:: +NonlinearSystemCompositional(Simulator& simulator, + const ModelParameters& param, + CompWellModel& wellModel, + const bool terminalOutput) + : ParentType(simulator, param, wellModel, terminalOutput) +{ + this->convergence_reports_.reserve(64); +} + +template +SimulatorReportSingle +NonlinearSystemCompositional:: +prepareStep(const SimulatorTimerInterface& timer) +{ + SimulatorReportSingle report; + Dune::Timer perfTimer; + perfTimer.start(); + + const int lastStepFailed = timer.lastStepFailed(); + if (this->grid_.comm().size() > 1 + && this->grid_.comm().max(lastStepFailed) != this->grid_.comm().min(lastStepFailed)) { + OPM_THROW(std::runtime_error, + "Misalignment of the parallel simulation run in prepareStep " + "- the previous step succeeded on some ranks but failed on others."); + } + + if (lastStepFailed) { + this->wellModel().restoreLastValidState(); + this->simulator_.model().updateFailed(); + } + else { + this->simulator_.model().advanceTimeLevel(); + } + + this->simulator_.setTime(timer.simulationTimeElapsed()); + this->simulator_.setTimeStepSize(timer.currentStepLength()); + + this->simulator_.problem().resetIterationForNewTimestep(); + this->simulator_.problem().beginTimeStep(); + + report.pre_post_time += perfTimer.stop(); + return report; +} + +template +void +NonlinearSystemCompositional:: +initialLinearization(SimulatorReportSingle& report, + const int minIter, + const int maxIter, + const SimulatorTimerInterface& timer) +{ + ParentType::initialLinearization(report, + minIter, + maxIter, + timer); + Dune::Timer perfTimer; + perfTimer.start(); + + auto convrep = ConvergenceReport{timer.simulationTimeElapsed()}; + const auto residualMetrics = this->reservoirResidualMetrics(); + const auto tolerance = this->simulator_.model().newtonMethod().tolerance(); + + for (int compIdx = 0; compIdx < numEq; ++compIdx) { + const std::array residual{residualMetrics[compIdx]}; + const std::array types{ + ConvergenceReport::ReservoirFailure::Type::MassBalance + }; + const std::array tolerances{tolerance}; + + this->addReservoirConvergenceMetrics( + convrep, + compIdx, + this->compNames_.name(compIdx), + residual, + types, + tolerances, + this->param_.max_residual_allowed_, + [this](const std::string& message) + { + if (this->terminal_output_) { + OpmLog::debug(message); + } + }); + } + + const auto severity = convrep.severityOfWorstFailure(); + const bool wellConverged = this->wellModel().getWellConvergence(); + report.converged = convrep.converged() && wellConverged && + this->simulator_.problem().iterationContext().iteration() >= minIter; + + this->convergence_reports_.back().report.push_back(std::move(convrep)); + report.update_time += perfTimer.stop(); + this->residual_norms_history_.push_back(residualMetrics); + + if (severity == ConvergenceReport::Severity::NotANumber) { + this->failureReport_ += report; + OPM_THROW_PROBLEM(NumericalProblem, "NaN residual found!"); + } + + if (severity == ConvergenceReport::Severity::TooLarge) { + this->failureReport_ += report; + OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); + } +} + +template +template +SimulatorReportSingle +NonlinearSystemCompositional:: +nonlinearIteration(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinearSolver) +{ + if (this->simulator_.problem().iterationContext().needsTimestepInit()) { + this->residual_norms_history_.clear(); + this->current_relaxation_ = 1.0; + this->dx_old_ = 0.0; + this->convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}}); + this->convergence_reports_.back().report.reserve(numEq); + } + + auto result = this->nonlinearIterationNewton(timer, nonlinearSolver); + this->simulator_.problem().advanceIteration(); + return result; +} + +template +template +SimulatorReportSingle +NonlinearSystemCompositional:: +nonlinearIterationNewton(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinearSolver) +{ + OPM_TIMEFUNCTION(); + + SimulatorReportSingle report; + Dune::Timer perfTimer; + + this->initialLinearization(report, + this->param_.newton_min_iter_, + this->param_.newton_max_iter_, + timer); + + if (!report.converged) { + perfTimer.reset(); + perfTimer.start(); + report.total_newton_iterations = 1; + + BVector x(this->simulator_.model().numGridDof()); + this->linear_solve_setup_time_ = 0.0; + + try { + auto& linearizer = this->simulator_.model().linearizer(); + linearizer.linearizeAuxiliaryEquations(); + linearizer.finalize(); + + this->solveJacobianSystem(x); + + report.linear_solve_setup_time += this->linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += this->linearIterationsLastSolve(); + } + catch (...) { + report.linear_solve_setup_time += this->linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += this->linearIterationsLastSolve(); + + this->failureReport_ += report; + throw; + } + + perfTimer.reset(); + perfTimer.start(); + + auto& model = this->simulator_.model(); + for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) { + model.auxiliaryModule(auxModIdx)->postSolve(x); + } + + if (this->param_.use_update_stabilization_) { + bool isOscillate = false; + bool isStagnate = false; + nonlinearSolver.detectOscillations(this->residual_norms_history_, + this->residual_norms_history_.size() - 1, + isOscillate, + isStagnate); + + if (isOscillate) { + this->current_relaxation_ -= nonlinearSolver.relaxIncrement(); + this->current_relaxation_ = std::max(this->current_relaxation_, nonlinearSolver.relaxMax()); + + if (this->terminalOutputEnabled()) { + OpmLog::info(" Oscillating behavior detected: Relaxation set to " + + std::to_string(this->current_relaxation_)); + } + } + + nonlinearSolver.stabilizeNonlinearUpdate(x, this->dx_old_, this->current_relaxation_); + } + + this->updateSolution(x); + report.update_time += perfTimer.stop(); + } + + return report; +} + +template +typename NonlinearSystemCompositional::Scalar +NonlinearSystemCompositional:: +relativeChange() const +{ + Scalar resultDelta = 0.0; + Scalar resultDenom = 0.0; + + const auto& elemMapper = this->simulator_.model().elementMapper(); + const auto& gridView = this->simulator_.gridView(); + + for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { + const unsigned globalElemIdx = elemMapper.index(elem); + const auto& priVarsNew = this->simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; + const auto& priVarsOld = this->simulator_.model().solution(/*timeIdx=*/1)[globalElemIdx]; + + for (int pvIdx = 0; pvIdx < static_cast(priVarsNew.size()); ++pvIdx) { + const auto delta = priVarsNew[pvIdx] - priVarsOld[pvIdx]; + resultDelta += delta * delta; + resultDenom += priVarsNew[pvIdx] * priVarsNew[pvIdx]; + } + } + + resultDelta = gridView.comm().sum(resultDelta); + resultDenom = gridView.comm().sum(resultDenom); + + return resultDenom > 0.0 ? resultDelta / resultDenom : 0.0; +} + +template +void +NonlinearSystemCompositional:: +solveJacobianSystem(BVector& x) +{ + auto& jacobian = this->simulator_.model().linearizer().jacobian(); + auto& residual = this->simulator_.model().linearizer().residual(); + auto& linSolver = this->simulator_.model().newtonMethod().linearSolver(); + + x = 0.0; + + Dune::Timer perfTimer; + perfTimer.start(); + linSolver.prepare(jacobian, residual); + this->linear_solve_setup_time_ = perfTimer.stop(); + linSolver.setResidual(residual); + linSolver.getResidual(residual); + linSolver.setMatrix(jacobian); + linSolver.solve(x); +} + +template +std::vector::Scalar> +NonlinearSystemCompositional:: +reservoirResidualMetrics() const +{ + const auto& model = this->simulator_.model(); + const auto& residual = model.linearizer().residual(); + const auto& constraintsMap = model.linearizer().constraintsMap(); + + std::vector residualMetrics(numEq, 0.0); + + for (unsigned dofIdx = 0; dofIdx < residual.size(); ++dofIdx) { + if (dofIdx >= model.numGridDof() || model.dofTotalVolume(dofIdx) <= 0.0) { + continue; + } + + if (constraintsMap.count(dofIdx) > 0) { + continue; + } + + const auto& localResidual = residual[dofIdx]; + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + residualMetrics[eqIdx] = std::max( + residualMetrics[eqIdx], + std::abs(localResidual[eqIdx] * model.eqWeight(dofIdx, eqIdx))); + } + } + + if (this->grid_.comm().size() > 1 && !residualMetrics.empty()) { + this->grid_.comm().max(residualMetrics.data(), residualMetrics.size()); + } + + return residualMetrics; +} + +} // namespace Opm + +#endif // OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/flow/BlackoilModelNldd.hpp b/opm/simulators/flow/NonlinearSystemNldd.hpp similarity index 96% rename from opm/simulators/flow/BlackoilModelNldd.hpp rename to opm/simulators/flow/NonlinearSystemNldd.hpp index 10ca0d0a2ad..e4f05b8a9c0 100644 --- a/opm/simulators/flow/BlackoilModelNldd.hpp +++ b/opm/simulators/flow/NonlinearSystemNldd.hpp @@ -21,8 +21,8 @@ along with OPM. If not, see . */ -#ifndef OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED -#define OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED +#ifndef OPM_NONLINEAR_SYSTEM_NLDD_HEADER_INCLUDED +#define OPM_NONLINEAR_SYSTEM_NLDD_HEADER_INCLUDED #include #include @@ -77,11 +77,11 @@ namespace Opm { -template class BlackoilModel; +template class NonlinearSystemBlackOilReservoir; -/// A NLDD implementation for three-phase black oil. +/// A NLDD implementation used by the reservoir nonlinear system. template -class BlackoilModelNldd +class NonlinearSystemNldd { public: using ElementContext = GetPropType; @@ -89,21 +89,19 @@ class BlackoilModelNldd using Grid = GetPropType; using Indices = GetPropType; using Scalar = GetPropType; - using ModelParameters = BlackoilModelParameters; + using ModelParameters = typename NonlinearSystemBlackOilReservoir::ModelParameters; using SolutionVector = GetPropType; - using BVector = typename BlackoilModel::BVector; + using BVector = typename NonlinearSystemBlackOilReservoir::BVector; using Domain = SubDomain; using ISTLSolverType = ISTLSolver; - using Mat = typename BlackoilModel::Mat; + using Mat = typename NonlinearSystemBlackOilReservoir::Mat; static constexpr int numEq = Indices::numEq; //! \brief The constructor sets up the subdomains. - //! \param model BlackOil model to solve for - //! \param param param Model parameters - //! \param compNames Names of the solution components - explicit BlackoilModelNldd(BlackoilModel& model) + //! \param model Owning nonlinear system to solve for + explicit NonlinearSystemNldd(NonlinearSystemBlackOilReservoir& model) : model_(model) , wellModel_(model.wellModel()) , rank_(model_.simulator().vanguard().grid().comm().rank()) @@ -849,22 +847,18 @@ class BlackoilModelNldd CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, CR::ReservoirFailure::Type::Cnv }; Scalar tol[2] = { tol_mb, tol_cnv }; - for (int ii : {0, 1}) { - if (std::isnan(res[ii])) { - report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); - logger.debug("NaN residual for " + model_.compNames().name(compIdx) + " equation."); - } else if (res[ii] > model_.param().max_residual_allowed_) { - report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); - logger.debug("Too large residual for " + model_.compNames().name(compIdx) + " equation."); - } else if (res[ii] < 0.0) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); - logger.debug("Negative residual for " + model_.compNames().name(compIdx) + " equation."); - } else if (res[ii] > tol[ii]) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); - } - - report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]); - } + model_.addReservoirConvergenceMetrics( + report, + compIdx, + model_.compNames().name(compIdx), + std::span{res}, + std::span{types}, + std::span{tol}, + model_.param().max_residual_allowed_, + [&logger](const std::string& message) + { + logger.debug(message); + }); } // Output of residuals. If converged at initial state, log nothing. @@ -1218,7 +1212,7 @@ class BlackoilModelNldd return false; } - BlackoilModel& model_; //!< Reference to model + NonlinearSystemBlackOilReservoir& model_; //!< Reference to model BlackoilWellModelNldd wellModel_; //!< NLDD well model adapter std::vector domains_; //!< Vector of subdomains std::vector> domain_matrices_; //!< Vector of matrix operator for each subdomain @@ -1235,4 +1229,4 @@ class BlackoilModelNldd } // namespace Opm -#endif // OPM_BLACKOILMODEL_NLDD_HEADER_INCLUDED +#endif // OPM_NONLINEAR_SYSTEM_NLDD_HEADER_INCLUDED diff --git a/opm/simulators/flow/NonlinearSystem_impl.hpp b/opm/simulators/flow/NonlinearSystem_impl.hpp new file mode 100644 index 00000000000..897d6731189 --- /dev/null +++ b/opm/simulators/flow/NonlinearSystem_impl.hpp @@ -0,0 +1,307 @@ +#ifndef OPM_FLOW_NONLINEAR_SYSTEM_IMPL_HEADER_INCLUDED +#define OPM_FLOW_NONLINEAR_SYSTEM_IMPL_HEADER_INCLUDED +/* + Copyright 2026, 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 . +*/ +#ifndef OPM_FLOW_NONLINEAR_SYSTEM_HEADER_INCLUDED +#include +#include +#endif + +#include + +#include +#include + +#include +#include +#include +#include + +namespace Opm { + +template +SimulatorReportSingle +NonlinearSystem:: +assembleReservoir(const SimulatorTimerInterface&) +{ + return assembleReservoir(well_model_); +} + +template +void +NonlinearSystem:: +updateTUNING(const Tuning& tuning) +{ + applyTUNING(param_, tuning); +} + +template +void +NonlinearSystem:: +updateTUNINGDP(const TuningDp& tuning_dp) +{ + applyTUNINGDP(param_, tuning_dp); +} + +template +void +NonlinearSystem:: +updateSolution(const GlobalEqVector& dx) +{ + OPM_TIMEBLOCK(updateSolution); + + const bool shouldStore = shouldStoreSolutionUpdate(); + if (shouldStore) { + prepareSolutionUpdate(); + } + + auto& newtonMethod = simulator_.model().newtonMethod(); + auto& solution = simulator_.model().solution(/*timeIdx=*/0); + + newtonMethod.applyUpdate(/*nextSolution=*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx); + + { + OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities); + simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + + if (shouldStore) { + storeSolutionUpdate(dx); + } +} + +template +template +void +NonlinearSystem:: +addReservoirConvergenceMetrics(ConvergenceReport& report, + const int componentIdx, + const std::string_view componentName, + const std::span residuals, + const std::span types, + const std::span tolerances, + const Scalar maxResidualAllowed, + LogFailure&& logFailure) const +{ + if (residuals.size() != types.size() || residuals.size() != tolerances.size()) { + OPM_THROW(std::logic_error, "Mismatched reservoir convergence metric sizes."); + } + + using CR = ConvergenceReport; + for (std::size_t metricIdx = 0; metricIdx < residuals.size(); ++metricIdx) { + const auto residual = residuals[metricIdx]; + const auto type = types[metricIdx]; + const auto tolerance = tolerances[metricIdx]; + + if (std::isnan(residual)) { + report.setReservoirFailed({type, CR::Severity::NotANumber, componentIdx}); + std::forward(logFailure)("NaN residual for " + std::string(componentName) + " equation."); + } + else if (residual > maxResidualAllowed) { + report.setReservoirFailed({type, CR::Severity::TooLarge, componentIdx}); + std::forward(logFailure)("Too large residual for " + std::string(componentName) + " equation."); + } + else if (residual < 0.0) { + report.setReservoirFailed({type, CR::Severity::Normal, componentIdx}); + std::forward(logFailure)("Negative residual for " + std::string(componentName) + " equation."); + } + else if (residual > tolerance) { + report.setReservoirFailed({type, CR::Severity::Normal, componentIdx}); + } + + report.setReservoirConvergenceMetric(type, componentIdx, residual, tolerance); + } +} + +template +NonlinearSystem:: +NonlinearSystem(Simulator& simulator, + const ModelParameters& param, + WellModel& wellModel, + const bool terminal_output) + : simulator_(simulator) + , grid_(simulator_.vanguard().grid()) + , terminal_output_(terminal_output) + , param_(param) + , well_model_(wellModel) + , current_relaxation_(1.0) + , dx_old_(simulator_.model().numGridDof()) +{} + +template +void +NonlinearSystem:: +initialLinearization(SimulatorReportSingle& report, + const int, + const int, + const SimulatorTimerInterface&) +{ + failureReport_ = SimulatorReportSingle(); + + Dune::Timer perfTimer; + perfTimer.start(); + report.total_linearizations = 1; + + try { + report += this->assembleReservoir(well_model_); + report.assemble_time += perfTimer.stop(); + + // Mark timestep initialized after assembling, because well models can use + // needsTimestepInit() to trigger per-step setup during assembly. + simulator_.problem().markTimestepInitialized(); + } + catch (...) { + report.assemble_time += perfTimer.stop(); + failureReport_ += report; + throw; + } +} + +template +SimulatorReportSingle +NonlinearSystem:: +prepareStep(const SimulatorTimerInterface& timer) +{ + SimulatorReportSingle report; + Dune::Timer perfTimer; + perfTimer.start(); + + const int lastStepFailed = timer.lastStepFailed(); + if (grid_.comm().size() > 1 && grid_.comm().max(lastStepFailed) != grid_.comm().min(lastStepFailed)) { + OPM_THROW(std::runtime_error, + "Misalignment of the parallel simulation run in prepareStep " + "- the previous step succeeded on some ranks but failed on others."); + } + + if (lastStepFailed) { + simulator_.model().updateFailed(); + } + else { + simulator_.model().advanceTimeLevel(); + } + + // The model still needs the report-step time context even though flow owns time stepping. + simulator_.setTime(timer.simulationTimeElapsed()); + simulator_.setTimeStepSize(timer.currentStepLength()); + + simulator_.problem().resetIterationForNewTimestep(); + simulator_.problem().beginTimeStep(); + + report.pre_post_time += perfTimer.stop(); + return report; +} + +template +template +SimulatorReportSingle +NonlinearSystem:: +assembleReservoir(WellModelType& wellModel) +{ + simulator_.problem().beginIteration(); + simulator_.model().linearizer().linearizeDomain(); + simulator_.problem().endIteration(); + return wellModel.lastReport(); +} + +template +template +void +NonlinearSystem:: +applyTUNING(ModelParametersType& param, + const Tuning& tuning) +{ + param.tolerance_cnv_ = tuning.TRGCNV; + param.tolerance_cnv_relaxed_ = tuning.XXXCNV; + param.tolerance_mb_ = tuning.TRGMBE; + param.tolerance_mb_relaxed_ = tuning.XXXMBE; + param.newton_max_iter_ = tuning.NEWTMX; + param.newton_min_iter_ = tuning.NEWTMN; +} + +template +template +void +NonlinearSystem:: +applyTUNINGDP(ModelParametersType& param, + const TuningDp& tuning_dp) +{ + param.tolerance_max_dp_ = tuning_dp.TRGDDP; + param.tolerance_max_ds_ = tuning_dp.TRGDDS; + param.tolerance_max_drs_ = tuning_dp.TRGDDRS; + param.tolerance_max_drv_ = tuning_dp.TRGDDRV; +} + +template +template +std::tuple +NonlinearSystem:: +convergenceReduction(Parallel::Communication comm, + const ValueType primaryVolumeLocal, + const ValueType secondaryVolumeLocal, + std::vector& sumValues, + std::vector& maxValues, + std::vector& averagedValues) +{ + OPM_TIMEBLOCK(convergenceReduction); + + ValueType primaryVolume = primaryVolumeLocal; + ValueType secondaryVolume = secondaryVolumeLocal; + + if (comm.size() > 1) { + std::vector sumBuffer; + std::vector maxBuffer; + const int numComp = averagedValues.size(); + sumBuffer.reserve(2 * numComp + 2); + maxBuffer.reserve(numComp); + + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + sumBuffer.push_back(averagedValues[compIdx]); + sumBuffer.push_back(sumValues[compIdx]); + maxBuffer.push_back(maxValues[compIdx]); + } + + sumBuffer.push_back(primaryVolume); + sumBuffer.push_back(secondaryVolume); + + comm.sum(sumBuffer.data(), sumBuffer.size()); + comm.max(maxBuffer.data(), maxBuffer.size()); + + for (int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx) { + averagedValues[compIdx] = sumBuffer[buffIdx]; + ++buffIdx; + sumValues[compIdx] = sumBuffer[buffIdx]; + } + + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + maxValues[compIdx] = maxBuffer[compIdx]; + } + + primaryVolume = sumBuffer[sumBuffer.size() - 2]; + secondaryVolume = sumBuffer.back(); + } + + return {primaryVolume, secondaryVolume}; +} + +} // namespace Opm + +#endif // OPM_FLOW_NONLINEAR_SYSTEM_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index a57a3b6d4df..e0d0a85eaed 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include @@ -53,6 +54,9 @@ #include #include +#include +#include +#include #include #include #include @@ -91,6 +95,12 @@ class OutputCompositionalModule : public GenericOutputBlackoilModule OutputCompositionalModule(const Simulator& simulator, const SummaryConfig& smryCfg, @@ -182,6 +192,92 @@ class OutputCompositionalModule : public GenericOutputBlackoilModule fipSched; + if (reportStepNum > 0) { + const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get(); + fipSched = std::make_unique(rpt); + } + + const FIPConfig& fipc = reportStepNum == 0 + ? this->eclState_.getEclipseConfig().fip() + : *fipSched; + + if (!substep && !this->forceDisableFipOutput_ && fipc.output(FIPConfig::OutputField::FIELD)) { + this->logOutput_.timeStamp("BALANCE", elapsed, reportStepNum, currentDate); + + const auto& initial_inplace = *this->initialInplace(); + this->logOutput_.fip(inplace, initial_inplace, ""); + + if (fipc.output(FIPConfig::OutputField::FIPNUM)) { + this->logOutput_.fip(inplace, initial_inplace, "FIPNUM"); + + if (fipc.output(FIPConfig::OutputField::RESV)) { + this->logOutput_.fipResv(inplace, "FIPNUM"); + } + } + + if (fipc.output(FIPConfig::OutputField::FIP)) { + for (const auto& reg : this->regions_) { + if (reg.first != "FIPNUM") { + std::ostringstream ss; + ss << "BAL" << reg.first.substr(3); + this->logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate); + this->logOutput_.fip(inplace, initial_inplace, reg.first); + + if (fipc.output(FIPConfig::OutputField::RESV)) { + this->logOutput_.fipResv(inplace, reg.first); + } + } + } + } + } + } + + void outputFipAndResvLogToCSV(const std::size_t reportStepNum, + const bool substep, + const Parallel::Communication& comm) + { + if (comm.rank() != 0) { + return; + } + + if ((reportStepNum == 0) && (!substep) && + (this->schedule_.initialReportConfiguration().has_value()) && + (this->schedule_.initialReportConfiguration()->contains("CSVFIP"))) { + + std::ostringstream csv_stream; + + this->logOutput_.csv_header(csv_stream); + + const auto& initial_inplace = *this->initialInplace(); + + this->logOutput_.fip_csv(csv_stream, initial_inplace, "FIPNUM"); + + for (const auto& reg : this->regions_) { + if (reg.first != "FIPNUM") { + this->logOutput_.fip_csv(csv_stream, initial_inplace, reg.first); + } + } + + const IOConfig& io = this->eclState_.getIOConfig(); + auto csv_fname = io.getOutputDir() + "/" + io.getBaseName() + ".CSV"; + + std::ofstream outputFile(csv_fname); + outputFile << csv_stream.str(); + } + } + //! \brief Setup list of active element-level data extractors void setupExtractors(const bool /*isSubStep*/, const std::size_t /*reportStepNum*/) diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp b/opm/simulators/flow/SimulatorFullyImplicit.cpp similarity index 98% rename from opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp rename to opm/simulators/flow/SimulatorFullyImplicit.cpp index 7223aa51774..f540693e6ec 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/simulators/flow/SimulatorFullyImplicit.cpp @@ -20,7 +20,7 @@ */ #include -#include +#include #include diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicit.hpp similarity index 96% rename from opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp rename to opm/simulators/flow/SimulatorFullyImplicit.hpp index e8e226fa930..281e171a6cd 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicit.hpp @@ -19,8 +19,8 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_HEADER_INCLUDED -#define OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_HEADER_INCLUDED +#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_HEADER_INCLUDED +#define OPM_SIMULATOR_FULLY_IMPLICIT_HEADER_INCLUDED #include #include @@ -37,7 +37,7 @@ #include #include -#include +#include #include #include #include @@ -83,7 +83,7 @@ void logTuning(const Tuning& tuning); namespace Opm { -/** \brief Top-level driver for a fully implicit black-oil simulation. +/** \brief Top-level driver for a fully implicit flow simulation. * * Owns the per-report-step loop: \ref run repeatedly invokes * \ref runStep until `timer.done()` is reached. Each \ref runStep @@ -110,7 +110,7 @@ namespace Opm { * program startup. */ template -class SimulatorFullyImplicitBlackoil : private SerializableSim +class SimulatorFullyImplicit : private SerializableSim { protected: struct MPI_Comm_Deleter; @@ -135,7 +135,7 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim using Solver = NonlinearSolver; using ModelParameters = typename Model::ModelParameters; using SolverParameters = typename Solver::SolverParameters; - using WellModel = BlackoilWellModel; + using WellModel = GetPropType; /** \brief Construct from the surrounding eWoms `Simulator`. * @@ -146,10 +146,10 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim * * \param simulator The surrounding eWoms simulator; observed, not owned. */ - explicit SimulatorFullyImplicitBlackoil(Simulator& simulator); + explicit SimulatorFullyImplicit(Simulator& simulator); /// Ends the convergence-output thread cleanly on all ranks. - ~SimulatorFullyImplicitBlackoil() override; + ~SimulatorFullyImplicit() override; /** \brief Register all parameters consumed by this class and its * major collaborators. @@ -372,6 +372,6 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim } // namespace Opm -#include +#include -#endif // OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_HEADER_INCLUDED +#endif // OPM_SIMULATOR_FULLY_IMPLICIT_HEADER_INCLUDED diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/simulators/flow/SimulatorFullyImplicit_impl.hpp similarity index 94% rename from opm/simulators/flow/SimulatorFullyImplicitBlackoil_impl.hpp rename to opm/simulators/flow/SimulatorFullyImplicit_impl.hpp index 99ac730d47b..0f3194f617d 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicit_impl.hpp @@ -19,13 +19,13 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_IMPL_HEADER_INCLUDED -#define OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_IMPL_HEADER_INCLUDED +#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED +#define OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED // Improve IDE experience -#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_HEADER_INCLUDED +#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_HEADER_INCLUDED #include -#include +#include #endif #include @@ -42,8 +42,8 @@ namespace Opm { template -SimulatorFullyImplicitBlackoil:: -SimulatorFullyImplicitBlackoil(Simulator& simulator) +SimulatorFullyImplicit:: +SimulatorFullyImplicit(Simulator& simulator) : simulator_(simulator) , serializer_(*this, FlowGenericVanguard::comm(), @@ -74,8 +74,8 @@ SimulatorFullyImplicitBlackoil(Simulator& simulator) } template -SimulatorFullyImplicitBlackoil:: -~SimulatorFullyImplicitBlackoil() +SimulatorFullyImplicit:: +~SimulatorFullyImplicit() { // Safe to call on all ranks, not just the I/O rank. convergence_output_.endThread(); @@ -83,7 +83,7 @@ SimulatorFullyImplicitBlackoil:: template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: registerParameters() { ModelParameters::registerParameters(); @@ -98,14 +98,14 @@ registerParameters() #ifdef RESERVOIR_COUPLING_ENABLED template SimulatorReport -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: run(SimulatorTimer& timer, int argc, char** argv) { init(timer, argc, argv); #else template SimulatorReport -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: run(SimulatorTimer& timer) { init(timer); @@ -142,7 +142,7 @@ run(SimulatorTimer& timer) #ifdef RESERVOIR_COUPLING_ENABLED template bool -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: checkRunningAsReservoirCouplingMaster() { for (std::size_t report_step = 0; report_step < this->schedule().size(); ++report_step) { @@ -170,7 +170,7 @@ checkRunningAsReservoirCouplingMaster() #ifdef RESERVOIR_COUPLING_ENABLED template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: init(const SimulatorTimer& timer, int argc, char** argv) { auto slave_mode = Parameters::Get(); @@ -200,7 +200,7 @@ init(const SimulatorTimer& timer, int argc, char** argv) #else template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: init(const SimulatorTimer& timer) { #endif @@ -236,7 +236,7 @@ init(const SimulatorTimer& timer) template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: updateTUNING(const Tuning& tuning) { modelParam_.tolerance_cnv_ = tuning.TRGCNV; @@ -252,7 +252,7 @@ updateTUNING(const Tuning& tuning) template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: updateTUNINGDP(const TuningDp& tuning_dp) { // NOTE: If TUNINGDP item is _not_ set it should be 0.0 @@ -275,7 +275,7 @@ updateTUNINGDP(const TuningDp& tuning_dp) template bool -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: runStep(SimulatorTimer& timer) { if (schedule().exitStatus().has_value()) { @@ -491,7 +491,7 @@ runStep(SimulatorTimer& timer) template SimulatorReport -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: finalize() { // make sure all output is written to disk before run is finished @@ -514,7 +514,7 @@ finalize() template template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: serializeOp(Serializer& serializer) { serializer(simulator_); @@ -524,7 +524,7 @@ serializeOp(Serializer& serializer) template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: loadState([[maybe_unused]] HDF5Serializer& serializer, [[maybe_unused]] const std::string& groupName) { @@ -535,7 +535,7 @@ loadState([[maybe_unused]] HDF5Serializer& serializer, template void -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: saveState([[maybe_unused]] HDF5Serializer& serializer, [[maybe_unused]] const std::string& groupName) const { @@ -546,7 +546,7 @@ saveState([[maybe_unused]] HDF5Serializer& serializer, template std::array -SimulatorFullyImplicitBlackoil:: +SimulatorFullyImplicit:: getHeader() const { std::ostringstream str; @@ -559,8 +559,8 @@ getHeader() const } template -std::unique_ptr::Solver> -SimulatorFullyImplicitBlackoil:: +std::unique_ptr::Solver> +SimulatorFullyImplicit:: createSolver(WellModel& wellModel) { auto model = std::make_unique(simulator_, @@ -591,4 +591,4 @@ createSolver(WellModel& wellModel) } // namespace Opm -#endif // OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_IMPL_HEADER_INCLUDED +#endif // OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/flow/TTagFlowProblemTPSA.hpp b/opm/simulators/flow/TTagFlowProblemTPSA.hpp index 45eed41f99e..41f8e628149 100644 --- a/opm/simulators/flow/TTagFlowProblemTPSA.hpp +++ b/opm/simulators/flow/TTagFlowProblemTPSA.hpp @@ -38,7 +38,7 @@ #include #include -#include +#include #include #include #include diff --git a/opm/simulators/flow/rescoup/ReservoirCouplingMaster.cpp b/opm/simulators/flow/rescoup/ReservoirCouplingMaster.cpp index 3bfc4d448b1..b87c33d2dd4 100644 --- a/opm/simulators/flow/rescoup/ReservoirCouplingMaster.cpp +++ b/opm/simulators/flow/rescoup/ReservoirCouplingMaster.cpp @@ -540,7 +540,7 @@ getMasterActivationDate_() const } } // NOTE: Consistency between SLAVES and GRUPMAST keywords has already been checked in - // init() in SimulatorFullyImplicitBlackoil.hpp + // init() in SimulatorFullyImplicit.hpp RCOUP_LOG_THROW(std::runtime_error, "Reservoir coupling: Failed to find master activation time: " "No SLAVES keyword found in schedule"); } diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index b6ec75edcee..b27c5a33a4d 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -95,7 +95,7 @@ class AdaptiveTimeStepping /** \brief Callback invoked at the start of each substep to apply * TUNING, NEXTSTEP (via ACTIONX), and WCYCLE updates. * - * Called once by `SimulatorFullyImplicitBlackoil::runStep` at the + * Called once by `SimulatorFullyImplicit::runStep` at the * start of a report step, and once per substep inside * \ref SubStepIteration::maybeUpdateTuningAndTimeStep_ (plus equivalent * sites in the reservoir-coupling loops). The callback may adjust @@ -280,7 +280,7 @@ class AdaptiveTimeStepping * \note The name is a historical artefact: by this point in the * report step, TUNING_CHANGE/TUNINGDP_CHANGE events have * already been cleared by the initial callback invocation - * in `SimulatorFullyImplicitBlackoil::runStep`, so what + * in `SimulatorFullyImplicit::runStep`, so what * actually runs here is the WCYCLE/NEXTSTEP branch of the * callback. */ @@ -350,7 +350,7 @@ class AdaptiveTimeStepping * * Splits the report step into a sequence of smaller substeps whose * lengths are chosen adaptively, invoking `solver` once per substep. - * Called from `SimulatorFullyImplicitBlackoil::runStep` in place of a + * Called from `SimulatorFullyImplicit::runStep` in place of a * direct `Solver::step()` call on the non-adaptive path. The actual * substep loop lives one layer down in \ref SubStepper::run(); this * method just constructs a \ref SubStepper and delegates. diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp index 024eb6cb26b..101b5d620d2 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -525,7 +525,7 @@ maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_ ); } -// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep() +// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicit::runStep() // It has to be called for each substep since TUNING might have been changed for next sub step due // to ACTIONX (via NEXTSTEP) or WCYCLE keywords. template @@ -1245,7 +1245,7 @@ maybeUpdateLastSubstepOfSyncTimestep_([[maybe_unused]] const double dt) #endif } -// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep() +// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicit::runStep() // It has to be called for each substep since TUNING might have been changed for next sub step due // to ACTIONX (via NEXTSTEP) or WCYCLE keywords. template @@ -1259,7 +1259,7 @@ maybeUpdateTuningAndTimeStep_() // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However, // the current definition of the maybeUpdateTuning_() callback is actually calling // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning - // see SimulatorFullyImplicitBlackoil::runStep() for more details. + // see SimulatorFullyImplicit::runStep() for more details. const auto old_value = suggestedNextTimestep_(); if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(), this->substep_timer_.currentStepLength(), diff --git a/opm/simulators/timestepping/ConvergenceReport.hpp b/opm/simulators/timestepping/ConvergenceReport.hpp index 563905738ab..d7de4f5d3bc 100644 --- a/opm/simulators/timestepping/ConvergenceReport.hpp +++ b/opm/simulators/timestepping/ConvergenceReport.hpp @@ -321,7 +321,7 @@ namespace Opm // Note regarding the CNV pore-volume split: We depend on the // fact that the quantities have already been aggregated across // all MPI ranks--see the implementation of member function - // BlackoilModel::getReservoirConvergence() for details--and are + // NonlinearSystemBlackOilReservoir::getReservoirConvergence() for details--and are // therefore equal on all ranks. Consequently, we simply assign // 'other's values here, if it is non-empty. Empty splits // typically come from well contributions. diff --git a/opm/simulators/utils/ComponentName.cpp b/opm/simulators/utils/ComponentName.cpp index 8ce0b2f5970..fbbe685e699 100644 --- a/opm/simulators/utils/ComponentName.cpp +++ b/opm/simulators/utils/ComponentName.cpp @@ -25,79 +25,51 @@ #include #include +#include #include #include #include +#include -#include +#include + +#include "ComponentName_impl.hpp" namespace Opm { -template -ComponentName::ComponentName() - : names_(Indices::numEq) -{ - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - names_[FluidSystem::canonicalToActiveCompIdx(canonicalCompIdx)] - = FluidSystem::componentName(canonicalCompIdx); - } - - if constexpr (Indices::enableSolvent) { - names_[Indices::solventSaturationIdx] = "Solvent"; - } - - if constexpr (Indices::enableExtbo) { - names_[Indices::zFractionIdx] = "ZFraction"; - } - - if constexpr (Indices::enablePolymer) { - names_[Indices::polymerConcentrationIdx] = "Polymer"; - } - - if constexpr (Indices::polymerMoleWeightIdx >= 0) { - assert(Indices::enablePolymer); - names_[Indices::polymerMoleWeightIdx] = "MolecularWeightP"; - } - - if constexpr (Indices::enableFullyImplicitThermal) { - names_[Indices::temperatureIdx] = "Energy"; - } - - if constexpr (Indices::numFoam == 1) { - names_[Indices::foamConcentrationIdx] = "Foam"; - } - - if constexpr (Indices::numBrine == 1) { - names_[Indices::saltConcentrationIdx] = "Brine"; - } - - if constexpr (Indices::enableMICP) { - names_[Indices::microbialConcentrationIdx] = "Microbes"; - names_[Indices::oxygenConcentrationIdx] = "Oxygen"; - names_[Indices::ureaConcentrationIdx] = "Urea"; - names_[Indices::biofilmVolumeFractionIdx] = "Biofilm"; - names_[Indices::calciteVolumeFractionIdx] = "Calcite"; - } - - if constexpr (Indices::enableBiofilm) { - names_[Indices::microbialConcentrationIdx] = "Microbes"; - names_[Indices::biofilmVolumeFractionIdx] = "Biofilm"; - } -} +#include +#define INSTANTIATE_COMPONENT_NAME_TYPE_TAG(T, N, W, TAG) \ +template class ComponentName, \ + FlashIndices, 0>>; -#include +#define INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, N, W) \ + INSTANTIATE_COMPONENT_NAME_TYPE_TAG(T, N, W, FlowCompProblem) + +#define INSTANTIATE_COMPONENT_NAME_COUNTS(T, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 2, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 3, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 4, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 5, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 6, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 7, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 8, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 9, W) \ + INSTANTIATE_COMPONENT_NAME_PROBLEMS(T, 10, W) INSTANTIATE_TYPE_INDICES(ComponentName, double) +INSTANTIATE_COMPONENT_NAME_COUNTS(double, false) +INSTANTIATE_COMPONENT_NAME_COUNTS(double, true) #if FLOW_INSTANTIATE_FLOAT INSTANTIATE_TYPE_INDICES(ComponentName, float) +INSTANTIATE_COMPONENT_NAME_COUNTS(float, false) +INSTANTIATE_COMPONENT_NAME_COUNTS(float, true) #endif +#undef INSTANTIATE_COMPONENT_NAME_COUNTS +#undef INSTANTIATE_COMPONENT_NAME_PROBLEMS +#undef INSTANTIATE_COMPONENT_NAME_TYPE_TAG + } // namespace Opm diff --git a/opm/simulators/utils/ComponentNameFlowExp.cpp b/opm/simulators/utils/ComponentNameFlowExp.cpp new file mode 100644 index 00000000000..fe0d46c86f4 --- /dev/null +++ b/opm/simulators/utils/ComponentNameFlowExp.cpp @@ -0,0 +1,41 @@ +#include + +#include + +#include + +#include + +#include + +#include "ComponentName_impl.hpp" + +namespace Opm { + +#define INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, N, W) \ +template class ComponentName, \ + FlashIndices, 0>>; + +#define INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS(T, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 2, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 3, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 4, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 5, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 6, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 7, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 8, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 9, W) \ + INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, 10, W) + +INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS(double, false) +INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS(double, true) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS(float, false) +INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS(float, true) +#endif + +#undef INSTANTIATE_FLOWEXP_COMPONENT_NAME_COUNTS +#undef INSTANTIATE_FLOWEXP_COMPONENT_NAME + +} // namespace Opm \ No newline at end of file diff --git a/opm/simulators/utils/ComponentName_impl.hpp b/opm/simulators/utils/ComponentName_impl.hpp new file mode 100644 index 00000000000..b40e650a73a --- /dev/null +++ b/opm/simulators/utils/ComponentName_impl.hpp @@ -0,0 +1,110 @@ +#ifndef OPM_COMPONENT_NAME_IMPL_HPP +#define OPM_COMPONENT_NAME_IMPL_HPP + +namespace Opm { + +template +ComponentName::ComponentName() + : names_(Indices::numEq) +{ + if constexpr (requires { + FluidSystem::solventComponentIndex(0u); + FluidSystem::canonicalToActiveCompIdx(0u); + }) { + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + names_[FluidSystem::canonicalToActiveCompIdx(canonicalCompIdx)] + = FluidSystem::componentName(canonicalCompIdx); + } + } + else { + constexpr auto numNamedComponents = + (Indices::numEq < FluidSystem::numComponents) ? Indices::numEq : FluidSystem::numComponents; + + for (unsigned compIdx = 0; compIdx < numNamedComponents; ++compIdx) { + names_[compIdx] = FluidSystem::componentName(compIdx); + } + + if constexpr (requires { Indices::waterEnabled; Indices::water0Idx; }) { + if constexpr (Indices::waterEnabled + && (Indices::water0Idx >= 0) + && (Indices::water0Idx < Indices::numEq)) { + names_[Indices::water0Idx] = "Water"; + } + } + } + + if constexpr (requires { Indices::enableSolvent; Indices::solventSaturationIdx; }) { + if constexpr (Indices::enableSolvent) { + names_[Indices::solventSaturationIdx] = "Solvent"; + } + } + + if constexpr (requires { Indices::enableExtbo; Indices::zFractionIdx; }) { + if constexpr (Indices::enableExtbo) { + names_[Indices::zFractionIdx] = "ZFraction"; + } + } + + if constexpr (requires { Indices::enablePolymer; Indices::polymerConcentrationIdx; }) { + if constexpr (Indices::enablePolymer) { + names_[Indices::polymerConcentrationIdx] = "Polymer"; + } + } + + if constexpr (requires { Indices::polymerMoleWeightIdx; }) { + if constexpr (Indices::polymerMoleWeightIdx >= 0) { + names_[Indices::polymerMoleWeightIdx] = "MolecularWeightP"; + } + } + + if constexpr (requires { Indices::enableFullyImplicitThermal; Indices::temperatureIdx; }) { + if constexpr (Indices::enableFullyImplicitThermal) { + names_[Indices::temperatureIdx] = "Energy"; + } + } + + if constexpr (requires { Indices::numFoam; Indices::foamConcentrationIdx; }) { + if constexpr (Indices::numFoam == 1) { + names_[Indices::foamConcentrationIdx] = "Foam"; + } + } + + if constexpr (requires { Indices::numBrine; Indices::saltConcentrationIdx; }) { + if constexpr (Indices::numBrine == 1) { + names_[Indices::saltConcentrationIdx] = "Brine"; + } + } + + if constexpr (requires { + Indices::enableMICP; + Indices::microbialConcentrationIdx; + Indices::oxygenConcentrationIdx; + Indices::ureaConcentrationIdx; + Indices::biofilmVolumeFractionIdx; + Indices::calciteVolumeFractionIdx; + }) { + if constexpr (Indices::enableMICP) { + names_[Indices::microbialConcentrationIdx] = "Microbes"; + names_[Indices::oxygenConcentrationIdx] = "Oxygen"; + names_[Indices::ureaConcentrationIdx] = "Urea"; + names_[Indices::biofilmVolumeFractionIdx] = "Biofilm"; + names_[Indices::calciteVolumeFractionIdx] = "Calcite"; + } + } + + if constexpr (requires { Indices::enableBiofilm; Indices::microbialConcentrationIdx; Indices::biofilmVolumeFractionIdx; }) { + if constexpr (Indices::enableBiofilm) { + names_[Indices::microbialConcentrationIdx] = "Microbes"; + names_[Indices::biofilmVolumeFractionIdx] = "Biofilm"; + } + } +} + +} // namespace Opm + +#endif \ No newline at end of file diff --git a/opm/simulators/wells/rescoup/RescoupProxy.hpp b/opm/simulators/wells/rescoup/RescoupProxy.hpp index 000796dcbfb..5d8aa7ca981 100644 --- a/opm/simulators/wells/rescoup/RescoupProxy.hpp +++ b/opm/simulators/wells/rescoup/RescoupProxy.hpp @@ -46,7 +46,7 @@ namespace ReservoirCoupling { /// providing mode queries that work regardless of MPI availability. This eliminates the need /// for #ifdef RESERVOIR_COUPLING_ENABLED guards in consumer code. /// -/// The proxy does not own the master/slave objects - they are owned by SimulatorFullyImplicitBlackoil. +/// The proxy does not own the master/slave objects - they are owned by SimulatorFullyImplicit. /// Pointers are set after construction via setMaster()/setSlave(). /// /// Master and slave modes are mutually exclusive - setting one clears the other. diff --git a/regressionTests.cmake b/regressionTests.cmake index dcc997edcce..947da790e3c 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -21,20 +21,14 @@ add_test_compareECLFiles( spe1 ) -add_test_compareECLFiles( - CASENAME - 1dcompositional - FILENAME - 1D_COMP - SIMULATOR - flowexp_comp - ABS_TOL - ${abs_tol} - REL_TOL - ${rel_tol} - DIR - compositional -) +add_test_compareECLFiles(CASENAME 1dcompositional + FILENAME 1D_COMP + SIMULATOR flowexp_comp + DEV_SIMULATOR flowexp_comp3_2p + REFERENCE_SIMULATOR flowexp_comp + ABS_TOL ${abs_tol} + REL_TOL ${rel_tol} + DIR compositional) add_test_compareECLFiles( CASENAME diff --git a/tests/rescoup/test_chopstep.cpp b/tests/rescoup/test_chopstep.cpp index 597b3a2d14e..e83989de19e 100644 --- a/tests/rescoup/test_chopstep.cpp +++ b/tests/rescoup/test_chopstep.cpp @@ -39,7 +39,7 @@ #include #include #include -#include +#include #include #if HAVE_DUNE_FEM diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index ee92bc531ba..9b7b02606d3 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -46,7 +46,7 @@ #include #include #include -#include +#include #include #include