OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-25 07:07:42 -05:00
Code Issues Packages Projects Releases Wiki Activity

Add standard-flow compositional support

Browse Source
This commit is contained in:
hnil committed 2026-06-11 20:08:03 +02:00
1 parent 95ab3ef446
commit f53f1c7198
30 files changed
+1702 -126
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists.txt
+58
View File
@@ -873,6 +873,64 @@ 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
$<TARGET_OBJECTS:moduleVersion>
LIBRARIES
opmsimulators
)
opm_add_executable(
TARGET
flowexp_comp3_2p
SOURCES
flowexperimental/comp/flowexp_comp3_2p_main.cpp
$<TARGET_OBJECTS:moduleVersion>
LIBRARIES
opmsimulators
)
opm_add_executable(
TARGET
flow_comp
ONLY_COMPILE
ALWAYS_ENABLE
DEPENDS
opmsimulators
LIBRARIES
opmsimulators
SOURCES
flowexperimental/comp/flow_comp.cpp
$<TARGET_OBJECTS:moduleVersion>
)
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
$<TARGET_OBJECTS:moduleVersion>
LIBRARIES
opmsimulators
)
opm_add_executable(
TARGET
flow_comp3_2p
SOURCES
flowexperimental/comp/flow_comp3_2p.cpp
$<TARGET_OBJECTS:moduleVersion>
LIBRARIES
opmsimulators
)
if(BUILD_FLOW_FLOAT_VARIANTS)
opm_add_executable(
TARGET
CMakeLists_files.cmake
+6
View File
@@ -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,6 +157,8 @@ 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
@@ -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
@@ -993,6 +997,8 @@ list (APPEND PUBLIC_HEADER_FILES
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
compareECLFiles.cmake
+21 -2
View File
@@ -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,7 +193,10 @@ function(add_test_compareSeparateECLFiles)
if(PARAM_IGNORE_EXTRA_KW)
list(APPEND DRIVER_ARGS -y ${PARAM_IGNORE_EXTRA_KW})
endif()
opm_add_test(${PARAM_PREFIX}_${PARAM_SIMULATOR}+${PARAM_CASENAME}
if(PARAM_REFERENCE_SIMULATOR)
list(APPEND DRIVER_ARGS -l ${PARAM_REFERENCE_SIMULATOR})
endif()
opm_add_test(${PARAM_PREFIX}_${PARAM_SIMULATOR}+${PARAM_CASENAME} NO_COMPILE
EXE_TARGET
${PARAM_SIMULATOR}
DRIVER_ARGS
@@ -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
flowexperimental/comp/flow_comp.cpp
+119
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/Parser/ParseContext.hpp>
#include <opm/input/eclipse/Parser/Parser.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
#include <fmt/format.h>
#include <array>
#include <cstdlib>
#include <tuple>
#include "flow_comp.hpp"
template <int compileTimeComponent>
std::tuple<bool, int>
runComponent(int runtimeComponent, bool water, int argc, char** argv)
{
if (runtimeComponent == compileTimeComponent) {
if (water)
return std::make_tuple(true, Opm::dispatchFlowComp<compileTimeComponent, true>(argc, argv));
else
return std::make_tuple(true, Opm::dispatchFlowComp<compileTimeComponent, false>(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 <int currentCompileTimeComponent, int nextComponent, int... components>
std::tuple<bool, int>
runComponent(int runtimecomponent, bool water, int argc, char** argv)
{
if (currentCompileTimeComponent == runtimecomponent) {
if (water)
return std::make_tuple(true, Opm::dispatchFlowComp<currentCompileTimeComponent, true>(argc, argv));
else
return std::make_tuple(true, Opm::dispatchFlowComp<currentCompileTimeComponent, false>(argc, argv));
}
return runComponent<nextComponent, components...>(runtimecomponent, water, argc, argv);
}
int
main(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::FlowCompProblem<0, true>;
Opm::registerEclTimeSteppingParameters<double>();
// 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_<TypeTag>(argc, const_cast<const char**>(argv), true, false, true, 0);
auto inputFilename
= Opm::FlowGenericVanguard::canonicalDeckPath(Opm::Parameters::Get<Opm::Parameters::EclDeckFileName>());
// 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<OPM_COMPILE_COMPONENTS_TEMPLATE_LIST>(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;
}
flowexperimental/comp/flow_comp.hpp
+301
View File
@@ -0,0 +1,301 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef FLOW_COMP_HPP
#define FLOW_COMP_HPP
#include <memory>
#include <opm/material/constraintsolvers/PTFlash.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/models/nonlinear/newtonmethod.hh>
#include <opm/models/ptflash/flashmodel.hh>
#include <opm/simulators/flow/Main.hpp>
#include <opm/simulators/flow/FlowProblemComp.hpp>
#include <opm/simulators/flow/FlowProblemCompProperties.hpp>
#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
#include <flowexperimental/comp/EmptyModel.hpp>
#include <flowexperimental/comp/wells/CompWellModel.hpp>
// // 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 numComp, bool EnableWater>
int dispatchFlowComp(int argc, char** argv);
} // namespace Opm
namespace Opm::Properties {
template <class TypeTag, class MyTypeTag>
struct DiscNewtonMethod;
template <class TypeTag, class MyTypeTag>
struct WellModel;
template <class TypeTag, class MyTypeTag>
struct TracerModel;
namespace TTag {
template<int NumComp, bool EnableWater>
struct FlowCompProblem {
using InheritsFrom = std::tuple<FlowBaseProblemComp, FlashModel>;
};
}
template<class TypeTag, int NumComp, bool EnableWater>
struct SparseMatrixAdapter<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
using Block = MatrixBlock<Scalar, numEq, numEq>;
public:
using type = typename Linear::IstlSparseMatrixAdapter<Block>;
};
// NOTE: Do NOT override NewtonMethod here.
// The inherited model configuration decides the solver implementation.
#if 0
template<class TypeTag>
struct SolidEnergyLaw<TypeTag, TTag::FlowCompProblem>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
public:
using EclThermalLawManager = ::Opm::EclThermalLawManager<Scalar, FluidSystem>;
using type = typename EclThermalLawManager::SolidEnergyLaw;
};
// Set the material law for thermal conduction
template<class TypeTag>
struct ThermalConductionLaw<TypeTag, TTag::FlowCompProblem>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
public:
using EclThermalLawManager = ::Opm::EclThermalLawManager<Scalar, FluidSystem>;
using type = typename EclThermalLawManager::ThermalConductionLaw;
};
template <class TypeTag>
struct SpatialDiscretizationSplice<TypeTag, TTag::FlowCompProblem>
{
using type = TTag::EcfvDiscretization;
};
template <class TypeTag>
struct LocalLinearizerSplice<TypeTag, TTag::FlowCompProblem>
{
using type = TTag::AutoDiffLocalLinearizer;
};
#endif
// Set the problem property
template <class TypeTag, int NumComp, bool EnableWater>
struct Problem<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>>
{
using type = FlowProblemComp<TypeTag>;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct AquiferModel<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
using type = EmptyModel<TypeTag>;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct WellModel<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
using type = CompWellModel<TypeTag>;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct TracerModel<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
using type = EmptyModel<TypeTag>;
};
template <class TypeTag, int NumComp, bool EnableWater>
struct FlashSolver<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
using type = Opm::PTFlash<Scalar, FluidSystem>;
};
template <class TypeTag, class MyTypeTag>
struct NumComp { using type = UndefinedProperty; };
// TODO: this is unfortunate, have to check why we need to hard-code it
template <class TypeTag, int NumComp_, bool EnableWater_>
struct NumComp<TypeTag, TTag::FlowCompProblem<NumComp_, EnableWater_>> {
static constexpr int value = NumComp_;
};
template <class TypeTag, class MyTypeTag>
struct EnableDummyWater { using type = UndefinedProperty; };
template <class TypeTag, int NumComp_, bool EnableWater_>
struct EnableDummyWater<TypeTag, TTag::FlowCompProblem<NumComp_, EnableWater_>> {
static constexpr bool value = EnableWater_;
};
#if 0
struct Temperature { using type = UndefinedProperty; };
template <class TypeTag>
struct Temperature<TypeTag, TTag::FlowCompProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 423.25;
};
#endif
template <class TypeTag, int NumComp_, bool EnableWater_>
struct FluidSystem<TypeTag, TTag::FlowCompProblem<NumComp_, EnableWater_>>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
static constexpr bool enable_water = getPropValue<TypeTag, Properties::EnableDummyWater>();
public:
using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp, enable_water>;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableMech<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableDisgasInWater<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct Stencil<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
public:
using type = EcfvStencil<Scalar, GridView>;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableApiTracking<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableSaltPrecipitation<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnablePolymerMW<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnablePolymer<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableDispersion<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableBrine<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableVapwat<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableSolvent<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableFoam<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableExtbo<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableBioeffects<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>> {
static constexpr bool value = false;
};
template<class TypeTag, int NumComp, bool EnableWater>
struct EnableDebuggingChecks<TypeTag, TTag::FlowCompProblem<NumComp, EnableWater>>{
static constexpr bool value = false;
};
// disable thermal flux boundaries by default
#if 0
template<class TypeTag>
struct EnableThermalFluxBoundaries<TypeTag, TTag::FlowCompProblem> {
static constexpr bool value = false;
};
#endif
} // namespace Opm::Properties
namespace Opm {
template<int numComp, bool EnableWater>
int dispatchFlowComp(int argc, char** argv)
{
using TypeTag = Properties::TTag::FlowCompProblem<numComp, EnableWater>;
auto mainObject = std::make_unique<Opm::Main>(argc, argv);
const auto ret = mainObject->runStatic<TypeTag>();
mainObject.reset();
return ret;
}
} // namespace Opm
#endif
flowexperimental/comp/flow_comp2_2p.cpp
+29
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
#include "flow_comp.hpp"
int main(int argc, char** argv)
{
return Opm::dispatchFlowComp<2, false>(argc, argv);
}
flowexperimental/comp/flow_comp3_2p.cpp
+29
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
#include "flow_comp.hpp"
int main(int argc, char** argv)
{
return Opm::dispatchFlowComp<3, false>(argc, argv);
}
flowexperimental/comp/flowexp_comp2_2p_main.cpp
+36
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
#include "flowexp_comp.hpp"
int main(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<2, false>;
Opm::registerEclTimeSteppingParameters<double>();
Opm::setupParameters_<TypeTag>(argc, const_cast<const char**>(argv), true, false, true, 0);
return Opm::start<TypeTag>(argc, argv, false);
}
flowexperimental/comp/flowexp_comp3_2p_main.cpp
+39
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/SimulatorFullyImplicit.hpp>
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
#include "flowexp_comp.hpp"
int main(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<3, false>;
Opm::Parameters::Register<Opm::Parameters::EnableAdaptiveTimeStepping>
("Use adaptive time stepping between report steps");
Opm::registerEclTimeSteppingParameters<double>();
Opm::setupParameters_<TypeTag>(argc, const_cast<const char**>(argv), true, false, true, 0);
return Opm::start<TypeTag>(argc, argv, false);
}
flowexperimental/comp/wells/CompWellModel.hpp
+27 -2
View File
@@ -23,6 +23,7 @@
#include <opm/output/data/Wells.hpp>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/input/eclipse/Schedule/Events.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
@@ -35,7 +36,10 @@
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <map>
#include <vector>
namespace Opm {
@@ -44,7 +48,7 @@ class Schedule;
struct NewtonIterationContext;
template<typename TypeTag>
class CompWellModel : WellConnectionAuxiliaryModule<TypeTag, CompWellModel<TypeTag>>
class CompWellModel : public WellConnectionAuxiliaryModule<TypeTag, CompWellModel<TypeTag>>
{
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@@ -96,10 +100,11 @@ public:
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 +131,19 @@ public:
auto begin() const { return well_container_.begin(); }
auto end() const { return well_container_.end(); }
const SimulatorReportSingle& lastReport() const { return last_report_; }
void prepareDeserialize(const int) {}
const std::map<std::string, double>& wellOpenTimes() const { return well_open_times_; }
const std::map<std::string, double>& 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 <class ReservoirCouplingSlave>
void setReservoirCouplingSlave(ReservoirCouplingSlave*) {}
template <class ReservoirCouplingMaster>
void setReservoirCouplingMaster(ReservoirCouplingMaster*) {}
bool getWellConvergence() const;
// the following functions are not used while added to avoid modifying WellConnectionAuxiliaryModule.hpp
@@ -146,6 +164,9 @@ private:
// we will need two to handle the changes between time stepping
CompWellState<FluidSystem> comp_well_states_;
// saved state at beginning of report step, used to restore on failed timestep
CompWellState<FluidSystem> last_valid_comp_well_states_;
// this is needed for parallel running, not all the wells will be in the same process
std::vector<Well> wells_ecl_;
std::vector<std::vector<CompConnectionData> > well_connection_data_;
@@ -155,6 +176,10 @@ private:
mutable BVector x_local_;
std::size_t local_num_cells_{0};
SimulatorReportSingle last_report_{};
std::map<std::string, double> well_open_times_;
std::map<std::string, double> well_close_times_;
WellGroupEvents report_step_start_events_{};
void createWellContainer();
void initWellContainer();
flowexperimental/comp/wells/CompWellModel_impl.hpp
+90 -3
View File
@@ -26,6 +26,8 @@
#include <flowexperimental/comp/wells/CompWellModel.hpp>
#endif
#include <algorithm>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
@@ -42,6 +44,7 @@ CompWellModel<TypeTag>::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<TypeTag>::
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 <typename TypeTag>
@@ -66,6 +95,24 @@ beginTimeStep()
initWellContainer();
}
template <typename TypeTag>
void
CompWellModel<TypeTag>::
restoreLastValidState()
{
comp_well_states_.copyDynamicStateFrom(last_valid_comp_well_states_);
}
template <typename TypeTag>
void
CompWellModel<TypeTag>::
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 <typename TypeTag>
void
CompWellModel<TypeTag>::
@@ -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<CompWell<TypeTag>>(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 <typename TypeTag>
bool
CompWellModel<TypeTag>::
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 <typename TypeTag>
bool
CompWellModel<TypeTag>::
flowexperimental/comp/wells/CompWellState.hpp
+9
View File
@@ -57,8 +57,17 @@ public:
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<SingleWellState> wells_;
flowexperimental/comp/wells/CompWellState_impl.hpp
+10 -3
View File
@@ -35,12 +35,19 @@ init(const std::vector<Well>& wells_ecl,
const std::vector<std::vector<Scalar>>& cell_mole_fractions,
const std::vector<std::vector<CompConnectionData> >& 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 <typename FluidSystem>
flowexperimental/comp/wells/CompWell_impl.hpp
+10
View File
@@ -17,6 +17,9 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <fmt/format.h>
#include <fmt/ranges.h>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
@@ -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;
}
flowexperimental/comp/wells/SingleCompWellState.hpp
+2
View File
@@ -94,6 +94,8 @@ public:
void update_injector_targets(const Well& well, const SummaryState& st);
void copyRuntimeStateFrom(const SingleCompWellState& other);
Scalar get_total_surface_rate() const;
};
flowexperimental/comp/wells/SingleCompWellState_impl.hpp
+16
View File
@@ -145,6 +145,22 @@ update_producer_targets(const Well& well,
}
}
template <typename FluidSystem>
void SingleCompWellState<FluidSystem>::
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 FluidSystem>
typename SingleCompWellState<FluidSystem>::Scalar
SingleCompWellState<FluidSystem>::
opm/models/utils/parametersystem.cpp
+17
View File
@@ -453,6 +453,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) {
opm/models/utils/parametersystem.hpp
+14
View File
@@ -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<Param>());
}
/*!
* \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 <class Param>
bool HideIfRegistered()
{
return detail::HideIfRegistered_(detail::getParamName<Param>());
}
/*!
* \brief Query whether parameter registration is open or not.
* \return True if registration is open, false otherwise
opm/simulators/flow/FlowGenericProblem.cpp
+55 -2
View File
@@ -29,6 +29,7 @@
#include <opm/grid/LookUpData.hh>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
@@ -40,29 +41,81 @@
namespace Opm {
using CpLeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>;
#define INSTANTIATE_TYPE(T) \
template class FlowGenericProblem< \
Dune::GridView< \
Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>, \
CpLeafGridView, \
BlackOilFluidSystem<T, BlackOilDefaultFluidSystemIndices>>;
#define INSTANTIATE_COMP_TYPE(T, N, W) \
template class FlowGenericProblem< \
CpLeafGridView, \
GenericOilGasWaterFluidSystem<T, N, W>>;
#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<Dune::CpGrid,
(Dune::PartitionIteratorType)4,
false>;
#define INSTANTIATE_DUNE_FEM_COMP_TYPE(T, N, W) \
template class FlowGenericProblem<GV, \
GenericOilGasWaterFluidSystem<T, N, W>>;
#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<GV,
BlackOilFluidSystem<double, BlackOilDefaultFluidSystemIndices>>;
INSTANTIATE_DUNE_FEM_FLOW_GENERIC_COMPONENTS(double)
#if FLOW_INSTANTIATE_FLOAT
template class FlowGenericProblem<GV,
BlackOilFluidSystem<float, BlackOilDefaultFluidSystemIndices>>;
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
opm/simulators/flow/FlowProblemComp.hpp
+12
View File
@@ -424,6 +424,18 @@ public:
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<class Serializer>
void serializeOp(Serializer& serializer)
opm/simulators/flow/FlowProblemCompProperties.hpp
+5
View File
@@ -33,6 +33,7 @@
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/flow/FlowBaseProblemProperties.hpp>
#include <opm/simulators/flow/NonlinearSystemCompositional.hpp>
#include <opm/material/thermal/EnergyModuleType.hpp>
#include <tuple>
@@ -51,6 +52,10 @@ struct FlowBaseProblemComp {
};
}
template<class TypeTag>
struct NonlinearSystem<TypeTag, TTag::FlowBaseProblemComp>
{ using type = NonlinearSystemCompositional<TypeTag>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::FlowBaseProblemComp>
opm/simulators/flow/FlowUtils.cpp
+44 -39
View File
@@ -215,47 +215,52 @@ void hideUnusedParameters()
// the default eWoms checkpoint/restart mechanism does not work with flow
Parameters::Hide<Parameters::RestartTime<Scalar>>();
Parameters::Hide<Parameters::RestartWritingInterval>();
// 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::EnableVtkOutput>())){
Parameters::Hide<Parameters::VtkWriteOilFormationVolumeFactor>();
Parameters::Hide<Parameters::VtkWriteOilSaturationPressure>();
Parameters::Hide<Parameters::VtkWriteOilVaporizationFactor>();
Parameters::Hide<Parameters::VtkWritePorosity>();
Parameters::Hide<Parameters::VtkWritePotentialGradients>();
Parameters::Hide<Parameters::VtkWritePressures>();
Parameters::Hide<Parameters::VtkWritePrimaryVars>();
Parameters::Hide<Parameters::VtkWritePrimaryVarsMeaning>();
Parameters::Hide<Parameters::VtkWriteProcessRank>();
Parameters::Hide<Parameters::VtkWriteRelativePermeabilities>();
Parameters::Hide<Parameters::VtkWriteSaturatedGasOilVaporizationFactor>();
Parameters::Hide<Parameters::VtkWriteSaturatedOilGasDissolutionFactor>();
Parameters::Hide<Parameters::VtkWriteSaturationRatios>();
Parameters::Hide<Parameters::VtkWriteSaturations>();
Parameters::Hide<Parameters::VtkWriteTemperature>();
Parameters::Hide<Parameters::VtkWriteViscosities>();
Parameters::Hide<Parameters::VtkWriteWaterFormationVolumeFactor>();
Parameters::Hide<Parameters::VtkWriteGasDissolutionFactor>();
Parameters::Hide<Parameters::VtkWriteGasFormationVolumeFactor>();
Parameters::Hide<Parameters::VtkWriteGasSaturationPressure>();
Parameters::Hide<Parameters::VtkWriteIntrinsicPermeabilities>();
Parameters::Hide<Parameters::VtkWriteTracerConcentration>();
Parameters::Hide<Parameters::VtkWriteExtrusionFactor>();
Parameters::Hide<Parameters::VtkWriteFilterVelocities>();
Parameters::Hide<Parameters::VtkWriteDensities>();
Parameters::Hide<Parameters::VtkWriteDofIndex>();
Parameters::Hide<Parameters::VtkWriteMobilities>();
//}
Parameters::Hide<Parameters::VtkWriteAverageMolarMasses>();
Parameters::Hide<Parameters::VtkWriteFugacities>();
Parameters::Hide<Parameters::VtkWriteFugacityCoeffs>();
Parameters::Hide<Parameters::VtkWriteMassFractions>();
Parameters::Hide<Parameters::VtkWriteMolarities>();
Parameters::Hide<Parameters::VtkWriteMoleFractions>();
Parameters::Hide<Parameters::VtkWriteTotalMassFractions>();
Parameters::Hide<Parameters::VtkWriteTotalMoleFractions>();
// 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::EnableVtkOutput>())){
Parameters::HideIfRegistered<Parameters::VtkWriteOilFormationVolumeFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteOilSaturationPressure>();
Parameters::HideIfRegistered<Parameters::VtkWriteOilVaporizationFactor>();
Parameters::HideIfRegistered<Parameters::VtkWritePorosity>();
Parameters::HideIfRegistered<Parameters::VtkWritePotentialGradients>();
Parameters::HideIfRegistered<Parameters::VtkWritePressures>();
Parameters::HideIfRegistered<Parameters::VtkWritePrimaryVars>();
Parameters::HideIfRegistered<Parameters::VtkWritePrimaryVarsMeaning>();
Parameters::HideIfRegistered<Parameters::VtkWriteProcessRank>();
Parameters::HideIfRegistered<Parameters::VtkWriteRelativePermeabilities>();
Parameters::HideIfRegistered<Parameters::VtkWriteSaturatedGasOilVaporizationFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteSaturatedOilGasDissolutionFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteSaturationRatios>();
Parameters::HideIfRegistered<Parameters::VtkWriteSaturations>();
Parameters::HideIfRegistered<Parameters::VtkWriteTemperature>();
Parameters::HideIfRegistered<Parameters::VtkWriteViscosities>();
Parameters::HideIfRegistered<Parameters::VtkWriteWaterFormationVolumeFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteGasDissolutionFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteGasFormationVolumeFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteGasSaturationPressure>();
Parameters::HideIfRegistered<Parameters::VtkWriteIntrinsicPermeabilities>();
Parameters::HideIfRegistered<Parameters::VtkWriteTracerConcentration>();
Parameters::HideIfRegistered<Parameters::VtkWriteExtrusionFactor>();
Parameters::HideIfRegistered<Parameters::VtkWriteFilterVelocities>();
Parameters::HideIfRegistered<Parameters::VtkWriteDensities>();
Parameters::HideIfRegistered<Parameters::VtkWriteDofIndex>();
Parameters::HideIfRegistered<Parameters::VtkWriteMobilities>();
//}
Parameters::HideIfRegistered<Parameters::VtkWriteAverageMolarMasses>();
Parameters::HideIfRegistered<Parameters::VtkWriteFugacities>();
Parameters::HideIfRegistered<Parameters::VtkWriteFugacityCoeffs>();
Parameters::HideIfRegistered<Parameters::VtkWriteMassFractions>();
Parameters::HideIfRegistered<Parameters::VtkWriteMolarities>();
Parameters::HideIfRegistered<Parameters::VtkWriteMoleFractions>();
Parameters::HideIfRegistered<Parameters::VtkWriteTotalMassFractions>();
Parameters::HideIfRegistered<Parameters::VtkWriteTotalMoleFractions>();
Parameters::HideIfRegistered<Parameters::VtkWriteTortuosities>();
Parameters::HideIfRegistered<Parameters::VtkWriteDiffusionCoefficients>();
Parameters::HideIfRegistered<Parameters::VtkWriteEffectiveDiffusionCoefficients>();
// hide average density option
Parameters::Hide<Parameters::UseAverageDensityMsWells>();
Parameters::HideIfRegistered<Parameters::UseAverageDensityMsWells>();
}
namespace {
opm/simulators/flow/NonlinearSolver.cpp
+6 -1
View File
@@ -120,7 +120,12 @@ using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Size>>;
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)
opm/simulators/flow/NonlinearSystemCompositional.hpp
+122
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED
#define OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <opm/simulators/flow/BlackoilModelParameters.hpp>
#include <opm/simulators/flow/NonlinearSystem.hpp>
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <flowexperimental/comp/wells/CompWellModel.hpp>
#include <filesystem>
#include <vector>
namespace Opm {
template <class TypeTag>
class NonlinearSystemCompositional : public NonlinearSystem<TypeTag>
{
public:
using ParentType = NonlinearSystem<TypeTag>;
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<TypeTag, Properties::SparseMatrixAdapter>;
using ModelParameters = BlackoilModelParameters<Scalar>;
static constexpr int numEq = Indices::numEq;
using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
NonlinearSystemCompositional(Simulator& simulator,
const ModelParameters& param,
CompWellModel<TypeTag>& wellModel,
bool terminalOutput);
SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer);
void initialLinearization(SimulatorReportSingle& report,
int minIter,
int maxIter,
const SimulatorTimerInterface& timer) override;
template <class NonlinearSolverType>
SimulatorReportSingle nonlinearIteration(const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinearSolver);
template <class NonlinearSolverType>
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<SimulatorReport>& domainAccumulatedReports() const
{
static const std::vector<SimulatorReport> emptyReports{};
return emptyReports;
}
void writeNonlinearIterationsPerCell(const std::filesystem::path&) const {}
template<class T>
std::vector<std::vector<Scalar>> computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
{ return computeFluidInPlace(fipnum); }
std::vector<std::vector<Scalar>> computeFluidInPlace(const std::vector<int>&) const
{ return {}; }
void writePartitions(const std::filesystem::path&) const {}
private:
std::vector<Scalar> reservoirResidualMetrics() const;
double linear_solve_setup_time_ = 0.0;
};
} // namespace Opm
#include <opm/simulators/flow/NonlinearSystemCompositional_impl.hpp>
#endif // OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED
opm/simulators/flow/NonlinearSystemCompositional_impl.hpp
+335
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_IMPL_HEADER_INCLUDED
#define OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_IMPL_HEADER_INCLUDED
#ifndef OPM_NONLINEAR_SYSTEM_COMPOSITIONAL_HEADER_INCLUDED
#include <config.h>
#include <opm/simulators/flow/NonlinearSystemCompositional.hpp>
#endif
#include <dune/common/timer.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <algorithm>
#include <array>
#include <cmath>
namespace Opm {
template <class TypeTag>
NonlinearSystemCompositional<TypeTag>::
NonlinearSystemCompositional(Simulator& simulator,
const ModelParameters& param,
CompWellModel<TypeTag>& wellModel,
const bool terminalOutput)
: ParentType(simulator, param, wellModel, terminalOutput)
{
this->convergence_reports_.reserve(64);
}
template <class TypeTag>
SimulatorReportSingle
NonlinearSystemCompositional<TypeTag>::
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 <class TypeTag>
void
NonlinearSystemCompositional<TypeTag>::
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<Scalar, 1> residual{residualMetrics[compIdx]};
const std::array<ConvergenceReport::ReservoirFailure::Type, 1> types{
ConvergenceReport::ReservoirFailure::Type::MassBalance
};
const std::array<Scalar, 1> 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 <class TypeTag>
template <class NonlinearSolverType>
SimulatorReportSingle
NonlinearSystemCompositional<TypeTag>::
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 <class TypeTag>
template <class NonlinearSolverType>
SimulatorReportSingle
NonlinearSystemCompositional<TypeTag>::
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 <class TypeTag>
typename NonlinearSystemCompositional<TypeTag>::Scalar
NonlinearSystemCompositional<TypeTag>::
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<int>(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 <class TypeTag>
void
NonlinearSystemCompositional<TypeTag>::
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 <class TypeTag>
std::vector<typename NonlinearSystemCompositional<TypeTag>::Scalar>
NonlinearSystemCompositional<TypeTag>::
reservoirResidualMetrics() const
{
const auto& model = this->simulator_.model();
const auto& residual = model.linearizer().residual();
const auto& constraintsMap = model.linearizer().constraintsMap();
std::vector<Scalar> 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
opm/simulators/flow/OutputCompositionalModule.hpp
+96
View File
@@ -42,6 +42,7 @@
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/io/vtkblackoilparams.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>
@@ -52,6 +53,9 @@
#include <algorithm>
#include <cstddef>
#include <fstream>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>
#include <type_traits>
@@ -90,6 +94,12 @@ class OutputCompositionalModule : public GenericOutputBlackoilModule<GetPropType
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
public:
static void registerParameters()
{
BaseType::registerParameters();
VtkBlackoilParams::registerParameters();
}
template <class CollectDataToIORankType>
OutputCompositionalModule(const Simulator& simulator,
const SummaryConfig& smryCfg,
@@ -181,6 +191,92 @@ public:
BaseType::assignToSolution(sol);
}
void outputFipAndResvLog(const Inplace& inplace,
const std::size_t reportStepNum,
double elapsed,
boost::posix_time::ptime currentDate,
const bool substep,
const Parallel::Communication& comm)
{
if (comm.rank() != 0) {
return;
}
std::unique_ptr<FIPConfig> fipSched;
if (reportStepNum > 0) {
const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
fipSched = std::make_unique<FIPConfig>(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*/)
opm/simulators/utils/ComponentName.cpp
+31 -60
View File
@@ -25,80 +25,51 @@
#include <opm/simulators/utils/ComponentName.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/models/blackoil/blackoilvariableandequationindices.hh>
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/models/ptflash/flashindices.hh>
#include <cassert>
#include <limits>
#include <flowexperimental/comp/flow_comp.hpp>
#include "ComponentName_impl.hpp"
namespace Opm {
template<class FluidSystem, class Indices>
ComponentName<FluidSystem,Indices>::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 != std::numeric_limits<unsigned>::max()) {
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 <opm/simulators/utils/InstantiationIndicesMacros.hpp>
#define INSTANTIATE_COMPONENT_NAME_TYPE_TAG(T, N, W, TAG) \
template class ComponentName<GenericOilGasWaterFluidSystem<T, N, W>, \
FlashIndices<Properties::TTag::TAG<N, W>, 0>>;
#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
opm/simulators/utils/ComponentNameFlowExp.cpp
+41
View File
@@ -0,0 +1,41 @@
#include <config.h>
#include <opm/simulators/utils/ComponentName.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/models/ptflash/flashindices.hh>
#include <flowexperimental/comp/flowexp_comp.hpp>
#include "ComponentName_impl.hpp"
namespace Opm {
#define INSTANTIATE_FLOWEXP_COMPONENT_NAME(T, N, W) \
template class ComponentName<GenericOilGasWaterFluidSystem<T, N, W>, \
FlashIndices<Properties::TTag::FlowExpCompProblem<N, W>, 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
opm/simulators/utils/ComponentName_impl.hpp
+114
View File
@@ -0,0 +1,114 @@
#ifndef OPM_COMPONENT_NAME_IMPL_HPP
#define OPM_COMPONENT_NAME_IMPL_HPP
#include <cassert>
#include <limits>
namespace Opm {
template<class FluidSystem, class Indices>
ComponentName<FluidSystem,Indices>::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 != std::numeric_limits<unsigned>::max()) {
assert(Indices::enablePolymer);
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
regressionTests.cmake
+8 -14
View File
@@ -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