diff --git a/CMakeLists.txt b/CMakeLists.txt index 2dc6e3e744..d9b34500b1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -665,6 +665,7 @@ macro(opm-simulators_targets_hook) target_sources(test_group_higher_constraints PRIVATE $) target_sources(test_RestartSerialization PRIVATE $) target_sources(test_glift1 PRIVATE $) + target_sources(test_tpsa_localresidual PRIVATE $) if(MPI_FOUND) target_sources(test_chopstep PRIVATE $) endif() @@ -674,6 +675,7 @@ macro(opm-simulators_targets_hook) blackoil_legacyassembly blackoil_nohyst blackoil_temp + blackoil_tpsa biofilm brine brine_energy @@ -689,6 +691,7 @@ macro(opm-simulators_targets_hook) gaswater_brine gaswater_dissolution gaswater_dissolution_diffuse + gaswater_dissolution_tpsa gaswater_energy gaswater_saltprec_energy gaswater_saltprec_vapwat @@ -700,6 +703,7 @@ macro(opm-simulators_targets_hook) oilwater_polymer_injectivity onephase onephase_energy + onephase_tpsa polymer solvent solvent_foam diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0f162f12ae..49efb05772 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -112,10 +112,12 @@ list (APPEND MAIN_SOURCE_FILES opm/models/io/vtkprimaryvarsparams.cpp opm/models/io/vtkptflashparams.cpp opm/models/io/vtktemperatureparams.cpp + opm/models/io/vtktpsaparams.cpp opm/models/io/restart.cpp opm/models/nonlinear/newtonmethodparams.cpp opm/models/parallel/tasklets.cpp opm/models/parallel/threadmanager.cpp + opm/models/tpsa/tpsanewtonmethodparams.cpp opm/models/utils/parametersystem.cpp opm/models/utils/simulatorutils.cpp opm/models/utils/terminal.cpp @@ -132,6 +134,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/EclGenericWriter.cpp opm/simulators/flow/ExtboContainer.cpp opm/simulators/flow/ExtraConvergenceOutputThread.cpp + opm/simulators/flow/FacePropertiesTPSA.cpp opm/simulators/flow/FIPContainer.cpp opm/simulators/flow/FlowGenericProblem.cpp opm/simulators/flow/FlowGenericVanguard.cpp @@ -189,6 +192,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/PreconditionerFactory7.cpp opm/simulators/linalg/PropertyTree.cpp opm/simulators/linalg/setupPropertyTree.cpp + opm/simulators/linalg/TPSALinearSolverParameters.cpp opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp opm/simulators/timestepping/AdaptiveTimeStepping.cpp opm/simulators/timestepping/ConvergenceReport.cpp @@ -520,6 +524,9 @@ list (APPEND TEST_SOURCE_FILES tests/test_stoppedwells.cpp tests/test_ThreePointHorizontalSatfuncConsistencyChecks.cpp tests/test_timer.cpp + tests/test_tpsa_face_properties.cpp + tests/test_tpsa_localresidual.cpp + tests/test_tpsa_primaryvariables.cpp tests/test_vfpproperties.cpp tests/test_WaterSatfuncConsistencyChecks.cpp tests/test_wellmodel.cpp @@ -630,6 +637,7 @@ list (APPEND TEST_DATA_FILES tests/equil_humidwetgas.DATA tests/equil_rsvd_and_rvvd.DATA tests/equil_rsvd_and_rvvd_and_rvwvd.DATA + tests/tpsa_ex.data tests/wetgas.DATA tests/satfuncEPS_B.DATA tests/wells_manager_data.data @@ -779,6 +787,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/models/discretization/common/linearizationtype.hh opm/models/discretization/common/restrictprolong.hh opm/models/discretization/common/tpfalinearizer.hh + opm/models/discretization/common/tpsalinearizer.hpp opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh opm/models/discretization/ecfv/ecfvdiscretization.hh opm/models/discretization/ecfv/ecfvgridcommhandlefactory.hh @@ -849,6 +858,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/models/io/vtktemperaturemodule.hpp opm/models/io/vtktemperatureparams.hpp opm/models/io/vtktensorfunction.hh + opm/models/io/vtktpsamodule.hpp + opm/models/io/vtktpsaparams.hpp opm/models/io/vtkvectorfunction.hh opm/models/ncp/ncpboundaryratevector.hh opm/models/ncp/ncpextensivequantities.hh @@ -896,6 +907,13 @@ list (APPEND PUBLIC_HEADER_FILES opm/models/richards/richardsprimaryvariables.hh opm/models/richards/richardsproperties.hh opm/models/richards/richardsratevector.hh + opm/models/tpsa/elasticityindices.hpp + opm/models/tpsa/elasticitylocalresidualtpsa.hpp + opm/models/tpsa/elasticityprimaryvariables.hpp + opm/models/tpsa/tpsabaseproperties.hpp + opm/models/tpsa/tpsamodel.hpp + opm/models/tpsa/tpsanewtonmethod.hpp + opm/models/tpsa/tpsanewtonmethodparams.hpp opm/models/utils/alignedallocator.hh opm/models/utils/basicparameters.hh opm/models/utils/basicproperties.hh @@ -925,6 +943,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/BlackoilModelNldd.hpp opm/simulators/flow/BlackoilModelParameters.hpp opm/simulators/flow/BlackoilModelProperties.hpp + opm/simulators/flow/BlackoilModelTPSA.hpp opm/simulators/flow/CO2H2Container.hpp opm/simulators/flow/CollectDataOnIORank.hpp opm/simulators/flow/CollectDataOnIORank_impl.hpp @@ -939,6 +958,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/EquilInitializer.hpp opm/simulators/flow/ExtboContainer.hpp opm/simulators/flow/ExtraConvergenceOutputThread.hpp + opm/simulators/flow/FacePropertiesTPSA.hpp + opm/simulators/flow/FacePropertiesTPSA_impl.hpp opm/simulators/flow/FemCpGridCompat.hpp opm/simulators/flow/FIBlackoilModel.hpp opm/simulators/flow/FIPContainer.hpp @@ -954,6 +975,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/FlowProblemComp.hpp opm/simulators/flow/FlowProblemCompProperties.hpp opm/simulators/flow/FlowProblemParameters.hpp + opm/simulators/flow/FlowProblemTPSA.hpp opm/simulators/flow/FlowsContainer.hpp opm/simulators/flow/FlowUtils.hpp opm/simulators/flow/FlowsData.hpp @@ -995,6 +1017,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/SolutionContainers.hpp opm/simulators/flow/SubDomain.hpp opm/simulators/flow/TTagFlowProblemTPFA.hpp + opm/simulators/flow/TTagFlowProblemTPSA.hpp opm/simulators/flow/TTagFlowProblemGasWater.hpp opm/simulators/flow/TTagFlowProblemOnePhase.hpp opm/simulators/flow/TracerContainer.hpp @@ -1047,6 +1070,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/is_gpu_operator.hpp opm/simulators/linalg/ISTLSolver.hpp opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp + opm/simulators/linalg/ISTLSolverTPSA.hpp opm/simulators/linalg/istlpreconditionerwrappers.hh opm/simulators/linalg/istlsolverwrappers.hh opm/simulators/linalg/istlsparsematrixadapter.hh @@ -1091,6 +1115,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/SmallDenseMatrixUtils.hpp opm/simulators/linalg/setupPropertyTree.hpp opm/simulators/linalg/superlubackend.hh + opm/simulators/linalg/TPSALinearSolverParameters.hpp opm/simulators/linalg/twolevelmethodcpr.hh opm/simulators/linalg/vertexborderlistfromgrid.hh opm/simulators/linalg/weightedresidreductioncriterion.hh diff --git a/flow/flow_blackoil_tpsa.cpp b/flow/flow_blackoil_tpsa.cpp new file mode 100644 index 0000000000..9f3a2015be --- /dev/null +++ b/flow/flow_blackoil_tpsa.cpp @@ -0,0 +1,104 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" + +#include +#include + +#include +#include +#include + +#include + + +namespace Opm::Properties { + +namespace TTag { + +struct FlowBlackOilProblemTPSA +{ + using InheritsFrom = std::tuple; +}; + +} // namespace Opm::Properties::TTag + +// /// +// Flow related properties +// /// +template +struct Linearizer +{ using type = TpfaLinearizer; }; + +template +struct LocalResidual +{ using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion +{ static constexpr bool value = false; }; + +template +struct AvoidElementContext +{ static constexpr bool value = true; }; + +// /// +// TPSA related properties +// /// +template +struct EnableMech +{ static constexpr bool value = true; }; + +template +struct Problem +{ using type = FlowProblemTPSA; }; + +template +struct NonlinearSystem +{ using type = BlackoilModelTPSA; }; + +} // namespace Opm::Properties + +namespace Opm { + +// ----------------- Main program ----------------- +int flowBlackoilTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + resetLocale(); + + FlowMain + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +int flowBlackoilTpsaMainStandalone(int argc, char** argv) +{ + using TypeTag = Properties::TTag::FlowBlackOilProblemTPSA; + auto mainObject = std::make_unique(argc, argv); + auto ret = mainObject->runStatic(); + // Destruct mainObject as the destructor calls MPI_Finalize! + mainObject.reset(); + return ret; +} + +} // namespace Opm diff --git a/flow/flow_blackoil_tpsa.hpp b/flow/flow_blackoil_tpsa.hpp new file mode 100644 index 0000000000..6198401362 --- /dev/null +++ b/flow/flow_blackoil_tpsa.hpp @@ -0,0 +1,35 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef FLOW_BLACKOIL_TPSA_HPP +#define FLOW_BLACKOIL_TPSA_HPP + + +namespace Opm { + +//! \brief Main function used in flow binary. +int flowBlackoilTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_blackoil_tpsa binary. +int flowBlackoilTpsaMainStandalone(int argc, char** argv); + +} // namespace Opm + +#endif diff --git a/flow/flow_blackoil_tpsa_main.cpp b/flow/flow_blackoil_tpsa_main.cpp new file mode 100644 index 0000000000..6800603678 --- /dev/null +++ b/flow/flow_blackoil_tpsa_main.cpp @@ -0,0 +1,28 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" +#include + + +int main(int argc, char** argv) +{ + return Opm::flowBlackoilTpsaMainStandalone(argc, argv); +} diff --git a/flow/flow_gaswater_dissolution_tpsa.cpp b/flow/flow_gaswater_dissolution_tpsa.cpp new file mode 100644 index 0000000000..c9d8094a26 --- /dev/null +++ b/flow/flow_gaswater_dissolution_tpsa.cpp @@ -0,0 +1,134 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" + +#include +#include +#include + +#include +#include +#include + +#include + + +namespace Opm::Properties { + +namespace TTag { + +struct FlowGasWaterDissolutionProblemTPSA +{ + using InheritsFrom = std::tuple; +}; + +} // namespace Opm::Properties::TTag + +// /// +// Flow related properties +// /// +template +struct Linearizer +{ using type = TpfaLinearizer; }; + +template +struct LocalResidual +{ using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion +{ static constexpr bool value = false; }; + +template +struct EnableDisgasInWater +{ static constexpr bool value = true; }; + +template +struct EnableVapwat +{ static constexpr bool value = true; }; + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::FlowProblem; + using FluidSystem = GetPropType; + static constexpr EnergyModules energyModuleType = getPropValue(); + static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal; + +public: + using type = BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + numEnergyVars, + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()>; +}; + +// /// +// TPSA related properties +// /// +template +struct EnableMech +{ static constexpr bool value = true; }; + +template +struct Problem +{ using type = FlowProblemTPSA; }; + +template +struct NonlinearSystem +{ using type = BlackoilModelTPSA; }; + +} // namespace Opm::Properties + +namespace Opm { + +// ----------------- Main program ----------------- +int flowGasWaterDissolutionTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + resetLocale(); + + FlowMain + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +int flowGasWaterDissolutionTpsaMainStandalone(int argc, char** argv) +{ + using TypeTag = Properties::TTag::FlowGasWaterDissolutionProblemTPSA; + auto mainObject = std::make_unique(argc, argv); + auto ret = mainObject->runStatic(); + // Destruct mainObject as the destructor calls MPI_Finalize! + mainObject.reset(); + return ret; +} + +} // namespace Opm diff --git a/flow/flow_gaswater_dissolution_tpsa.hpp b/flow/flow_gaswater_dissolution_tpsa.hpp new file mode 100644 index 0000000000..8db759716a --- /dev/null +++ b/flow/flow_gaswater_dissolution_tpsa.hpp @@ -0,0 +1,34 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef FLOW_GAS_WATER_DISSOLUTION_TPSA_HPP +#define FLOW_GAS_WATER_DISSOLUTION_TPSA_HPP + +namespace Opm { + +//! \brief Main function used in flow binary. +int flowGasWaterDissolutionTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_gaswater_dissolution_tpsa binary. +int flowGasWaterDissolutionTpsaMainStandalone(int argc, char** argv); + +} // namespace Opm + +#endif diff --git a/flow/flow_gaswater_dissolution_tpsa_main.cpp b/flow/flow_gaswater_dissolution_tpsa_main.cpp new file mode 100644 index 0000000000..af7663d1d5 --- /dev/null +++ b/flow/flow_gaswater_dissolution_tpsa_main.cpp @@ -0,0 +1,28 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" +#include + + +int main(int argc, char** argv) +{ + return Opm::flowGasWaterDissolutionTpsaMainStandalone(argc, argv); +} diff --git a/flow/flow_onephase_tpsa.cpp b/flow/flow_onephase_tpsa.cpp new file mode 100644 index 0000000000..deae40ba7a --- /dev/null +++ b/flow/flow_onephase_tpsa.cpp @@ -0,0 +1,125 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" + +#include +#include +#include + +#include +#include +#include + +#include + + +namespace Opm::Properties { + +namespace TTag { + +struct FlowWaterOnlyProblemTPSA +{ + using InheritsFrom = std::tuple; +}; + +} // namespace Opm::Properties::TTag + +// /// +// Flow related properties +// /// +template +struct Linearizer +{ using type = TpfaLinearizer; }; + +template +struct LocalResidual +{ using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion +{ static constexpr bool value = false; }; + +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::FlowProblem; + using FluidSystem = GetPropType; + static constexpr EnergyModules energyModuleType = getPropValue(); + static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal; + +public: + using type = BlackOilOnePhaseIndices(), + getPropValue(), + getPropValue(), + numEnergyVars, + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*enabledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()>; +}; + +// /// +// TPSA related properties +// /// +template +struct EnableMech +{ static constexpr bool value = true; }; + +template +struct Problem +{ using type = FlowProblemTPSA; }; + +template +struct NonlinearSystem +{ using type = BlackoilModelTPSA; }; + +} // namespace Opm::Properties + +namespace Opm { + +// ----------------- Main program ----------------- +int flowWaterOnlyTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + resetLocale(); + + FlowMain + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +int flowWaterOnlyTpsaMainStandalone(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::FlowWaterOnlyProblemTPSA; + auto mainObject = std::make_unique(argc, argv); + auto ret = mainObject->runStatic(); + // Destruct mainObject as the destructor calls MPI_Finalize! + mainObject.reset(); + return ret; +} + +} // namespace Opm diff --git a/flow/flow_onephase_tpsa.hpp b/flow/flow_onephase_tpsa.hpp new file mode 100644 index 0000000000..31876f736a --- /dev/null +++ b/flow/flow_onephase_tpsa.hpp @@ -0,0 +1,34 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef FLOW_ONEPHASE_TPSA_HPP +#define FLOW_ONEPHASE_TPSA_HPP + +namespace Opm { + +//! \brief Main functon used in main flow binary. +int flowWaterOnlyTpsaMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_onephase_tpsa binary. +int flowWaterOnlyTpsaMainStandalone(int argc, char** argv); + +} // namespace Opm + +#endif diff --git a/flow/flow_onephase_tpsa_main.cpp b/flow/flow_onephase_tpsa_main.cpp new file mode 100644 index 0000000000..ad8f3a55f6 --- /dev/null +++ b/flow/flow_onephase_tpsa_main.cpp @@ -0,0 +1,28 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include "config.h" +#include + + +int main(int argc, char** argv) +{ + return Opm::flowWaterOnlyTpsaMainStandalone(argc, argv); +} diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index e8e236e8eb..5fe2739ad8 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -108,6 +108,7 @@ class BlackOilIntensiveQuantities enum { enableDispersion = getPropValue() }; enum { enableConvectiveMixing = getPropValue() }; enum { enableBioeffects = getPropValue() }; + enum { enableMech = getPropValue() }; enum { enableMICP = Indices::enableMICP }; enum { numPhases = getPropValue() }; enum { waterCompIdx = FluidSystem::waterCompIdx }; @@ -588,6 +589,31 @@ public: const Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx); porosity_ *= (1.0 - Sp); } + + + // Geomechanical updates to porosity/pore volume + if constexpr (enableMech) { + // TPSA calculations + if (problem.simulator().vanguard().eclState().runspec().mechSolver().tpsa()) { + // TPSA compressibility term + const Scalar rockBiot = problem.rockBiotComp(globalSpaceIdx); + if (rockBiot > 0.0) { + const Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx); + Evaluation active_pressure; + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + active_pressure = fluidState_.pressure(oilPhaseIdx) - rockRefPressure; + } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){ + active_pressure = fluidState_.pressure(waterPhaseIdx) - rockRefPressure; + } else { + active_pressure = fluidState_.pressure(gasPhaseIdx) - rockRefPressure; + } + porosity_ += rockBiot * active_pressure; + } + + // TPSA coupling term, pore volume changes due to mechanics + porosity_ += problem.rockMechPoroChange(globalSpaceIdx, /*timeIdx=*/timeIdx); + } + } } void assertFiniteMembers() diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index c29ae6b7b2..254a628f0d 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/opm/models/blackoil/blackoilmodel.hh @@ -220,6 +220,10 @@ template struct EnableConvectiveMixing { static constexpr bool value = false; }; +template +struct EnableMech +{ static constexpr bool value = false; }; + template struct RunAssemblyOnGpu { static constexpr bool value = false; }; diff --git a/opm/models/blackoil/blackoilproblem.hh b/opm/models/blackoil/blackoilproblem.hh index c6c7ce0838..87a0ff488a 100644 --- a/opm/models/blackoil/blackoilproblem.hh +++ b/opm/models/blackoil/blackoilproblem.hh @@ -139,6 +139,66 @@ public: Scalar rockCompressibility(unsigned) const { return 0.0; } + /*! + * \brief Returns the porosity (i.e., pore volume) change due to geomechanics + */ + template + Scalar rockMechPoroChange(const Context&, + unsigned, + unsigned) const + { return 0.0; } + + /*! + * \brief Returns the porosity (i.e., pore volume) change due to geomechanics + */ + Scalar rockMechPoroChange(unsigned) const + { return 0.0; } + + /*! + * \brief Returns the additional compressibility of a cell due to poroelasticity + */ + template + Scalar rockBiotComp(const Context&, + unsigned, + unsigned) const + { return 0.0; } + + /*! + * \brief Returns the additional compressibility of a cell due to poroelasticity + */ + Scalar rockBiotComp(unsigned) const + { return 0.0; } + + /*! + * \brief Returns Lame's first parameter of a cell + */ + template + Scalar lame(const Context&, + unsigned, + unsigned) const + { return 0.0; } + + /*! + * \brief Returns Lame's first parameter of a cell + */ + Scalar lame(unsigned) const + { return 0.0; } + + /*! + * \brief Returns Biot coefficient of a cell + */ + template + Scalar biotCoeff(const Context&, + unsigned, + unsigned) const + { return 0.0; } + + /*! + * \brief Returns Biot coefficient of a cell + */ + Scalar biotCoeff(unsigned) const + { return 0.0; } + /*! * \brief Returns the reference pressure for rock the compressibility of a cell */ diff --git a/opm/models/discretization/common/tpsalinearizer.hpp b/opm/models/discretization/common/tpsalinearizer.hpp new file mode 100644 index 0000000000..f6d10f1011 --- /dev/null +++ b/opm/models/discretization/common/tpsalinearizer.hpp @@ -0,0 +1,745 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_LINEARIZER_HPP +#define TPSA_LINEARIZER_HPP + +#include + +#include + +#include + +#include + +#include + +#include +#include + +#include +#include +#include + + +namespace Opm { + +/*! +* \brief Linearizes TPSA equations and generates system matrix and residual for linear solver +*/ +template +class TpsaLinearizer +{ + using Constraints = GetPropType; // TODO: make TPSA constraints + using Evaluation = GetPropType; + using FlowModel = GetPropType; + using GeomechModel = GetPropType; + using GlobalEqVector = GetPropType; + using GridView = GetPropType; + using LocalResidual = GetPropType; + using Problem = GetPropType; + using Scalar = GetPropType; + using Simulator = GetPropType; + using SparseMatrixAdapter = GetPropType; + using Stencil = GetPropType; + + using MaterialState = MaterialStateTPSA; + + enum { numEq = getPropValue() }; + + using ADVectorBlock = Dune::FieldVector; + using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock; + using VectorBlock = Dune::FieldVector; + +public: + // /// + // Public functions + // /// + /*! + * \brief Constructor + */ + TpsaLinearizer() + { + simulatorPtr_ = nullptr; + } + + /*! + * \brief Initialize the linearizer + * + * \param simulator Simulator object + * + * At this point we can assume that all objects in the simulator have been allocated. We cannot assume that they are + * fully initialized, though. + */ + void init(Simulator& simulator) + { + simulatorPtr_ = &simulator; + eraseMatrix(); + } + + /*! + * \brief Causes the Jacobian matrix to be recreated from scratch before the next iteration. + * + * This method is usually called if the sparsity pattern has changed for some reason. (e.g. by modifications of the + * grid or changes of the auxiliary equations.) + */ + void eraseMatrix() + { + jacobian_.reset(); + } + + /*! + * \brief Finalize creation of Jacobian matrix and make ready for linear solver + */ + void finalize() + { + jacobian_->finalize(); + } + + /*! + * \brief Initializing and/or reset residual and Jacobian + * + * \param domain (Sub-)domain to reset + */ + template + void resetSystem_(const SubDomainType& domain) + { + if (!jacobian_) { + initFirstIteration_(); + } + for (int globI : domain.cells) { + residual_[globI] = 0.0; + jacobian_->clearRow(globI, 0.0); + } + } + + /*! + * \brief Linearize the non-linear system + */ + void linearize() + { + linearizeDomain(); + linearizeAuxiliaryEquations(); + } + + /*! + * \brief Linearize the non-linear system for the spatial domain + * + * The global Jacobian of the residual is assembled and the residual is evaluated for the current solution. The + * current state of affairs (esp. the previous and the current solutions) is represented by the model object. + */ + void linearizeDomain() + { + int succeeded; + try { + linearizeDomain(fullDomain_); + succeeded = 1; + } + catch (const std::exception& e) { + std::cout << "rank " << simulator_().gridView().comm().rank() + << " caught an exception while linearizing TPSA system:" << e.what() + << "\n" << std::flush; + succeeded = 0; + } + catch (...) { + std::cout << "rank " << simulator_().gridView().comm().rank() + << " caught an exception while linearizing TPSA system" + << "\n" << std::flush; + succeeded = 0; + } + succeeded = simulator_().gridView().comm().min(succeeded); + + if (!succeeded) { + throw NumericalProblem("A process did not succeed in linearizing the TPSA system"); + } + } + + /*! + * \brief Linearize the non-linear system for the spatial domain + * + * \param domain (Sub-)domain to linearize + * + * The global Jacobian of the residual is assembled and the residual is evaluated for the current solution. The + * current state of affairs (esp. the previous and the current solutions) is represented by the model object. + */ + template + void linearizeDomain(const SubDomainType& domain) + { + OPM_TIMEBLOCK(linearizeDomain); + + // We defer the initialization of the Jacobian matrix until here because the auxiliary modules usually assume + // the problem, model and grid to be fully initialized' + if (!jacobian_) { + initFirstIteration_(); + } + + if (domain.cells.size() == flowModel_().numTotalDof()) { + // We are on the full domain. + resetSystem_(); + } + else { + resetSystem_(domain); + } + + linearize_(domain); + } + + /*! + * \brief Linearize auxillary equation + */ + void linearizeAuxiliaryEquations() + { } + + /*! + * \brief Extract local residuals and sub-block Jacobians from locally computed AD residual + * + * \param res Residual to set + * \param bMat Block matrix from Jacobian to set + * \param resid Computed AD residual + */ + void setResAndJacobi(VectorBlock& res, MatrixBlock& bMat, const ADVectorBlock& resid) const + { + // Scalar local residual + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + res[eqIdx] = resid[eqIdx].value(); + } + + // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the residual function 'eqIdx' for the + // degree of freedom 'dofIdx' with regard to the focus variable 'pvIdx' of the degree of freedom 'focusDofIdx' + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) { + bMat[eqIdx][pvIdx] = resid[eqIdx].derivative(pvIdx); + } + } + } + + /*! + * \brief Update boundary condition information + */ + void updateBoundaryConditionData() + { + for (auto& bdyInfo : boundaryInfo_) { + // Get boundary information from problem + const auto [type, displacementAD] = problem_().mechBoundaryCondition(bdyInfo.cell, bdyInfo.dir); + + // Strip the unnecessary (and zero anyway) derivatives off displacement + std::vector displacement(3, 0.0); + for (std::size_t ii = 0; ii < displacement.size(); ++ii) { + displacement[ii] = displacementAD[ii].value(); + } + + // Update boundary information + bdyInfo.bcdata.type = type; + bdyInfo.bcdata.displacement = displacement; + } + } + + // /// + // Public get and set functions + // /// + /*! + * \brief Get Jacobian matrix + * + * \returns Reference to (sparse) Jacobian matrix object + */ + const SparseMatrixAdapter& jacobian() const + { return *jacobian_; } + + /*! + * \brief Get Jacobian matrix + * + * \returns Reference to (sparse) Jacobian matrix object + */ + SparseMatrixAdapter& jacobian() + { return *jacobian_; } + + /*! + * \brief Get residual vector + * + * \returns Reference to the residual + */ + const GlobalEqVector& residual() const + { return residual_; } + + /*! + * \brief Get residual vector + * + * \returns Reference to the residual + */ + GlobalEqVector& residual() + { return residual_; } + + /*! + * \brief Get linearization type + * + * \returns Reference to the linearization object + * + * The LinearizationType controls the scheme used and the focus time index. The default is fully implicit scheme, and + * focus index equal to 0, i.e. current time (end of step). + */ + const LinearizationType& getLinearizationType() const + { return linearizationType_; } + + /*! + * \brief Get constraints map + * + * \returns Constraints map + * + * \note No constraints implemented in TPSA + */ + std::map constraintsMap() const + { return {}; } + + /*! + * \brief Set linearization type + * + * \param linearizationType Linearization object + */ + void setLinearizationType(LinearizationType linearizationType) + { linearizationType_ = linearizationType; } + +private: + // /// + // Private functions + // /// + /*! + * \brief Construct the Jacobian matrix and create cell-neighbor information + * + * Sets up cell-neigbor and boundary information structs for easy calculation in linearize_(). In addition, sparsity + * patterns for Jacobian matrix is set up, + */ + void createMatrix_() + { + OPM_TIMEBLOCK(createMatrixTPSA); + + // If the Jacobian has been initialize before, we jump out + if (!neighborInfo_.empty()) { + return; + } + + // Init. the stencil + const auto& flowModel = flowModel_(); + Stencil stencil(gridView_(), flowModel.dofMapper()); + + // Build up sparsity patterns and neighboring information for Jacobian and linearization + std::vector> sparsityPattern(flowModel.numTotalDof()); + unsigned numCells = flowModel.numTotalDof(); + neighborInfo_.reserve(numCells, 6 * numCells); + std::vector loc_nbinfo; + for (const auto& elem : elements(gridView_())) { + // Loop over primary dofs in the element + stencil.update(elem); + + for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) { + // Build up neighboring information for curret primary dof + const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx); + loc_nbinfo.resize(stencil.numDof() - 1); + + for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { + // NOTE: NeighborInfo could/should be expanded with cell face parameters located in problem_() + // needed when computing face terms in LocalResidual + const unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); + sparsityPattern[myIdx].insert(neighborIdx); + if (dofIdx > 0) { + const auto scvfIdx = dofIdx - 1; + const auto& scvf = stencil.interiorFace(scvfIdx); + const Scalar area = scvf.area(); + loc_nbinfo[dofIdx - 1] = NeighborInfo{ neighborIdx, area, nullptr }; + } + } + + // Insert local neighbor info + neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); + + // Boundary condition information + unsigned bfIndex = 0; + for (const auto& intersection : intersections(gridView_(), elem)) { + if (intersection.boundary()) { + // Get boundary face direction + const auto& bf = stencil.boundaryFace(bfIndex); + const int dir_id = bf.dirId(); + + // Skip NNCs + if (dir_id < 0) { + continue; + } + + // Get boundary information from problem() + const auto [type, displacementAD] = problem_().mechBoundaryCondition(myIdx, dir_id); + + // Strip the unnecessary (and zero anyway) derivatives off displacement + std::vector displacement(3, 0.0); + for (std::size_t ii = 0; ii < displacement.size(); ++ii) { + displacement[ii] = displacementAD[ii].value(); + } + + // Insert boundary condition data in container + BoundaryConditionData bcdata { type, displacement, bfIndex, bf.area() }; + boundaryInfo_.push_back( { myIdx, dir_id, bfIndex, bcdata } ); + ++bfIndex; + continue; + } + if (!intersection.neighbor()) { + ++bfIndex; + continue; + } + } + } + } + + // Allocate Jacobian matrix and pointers to its sub-blocks + jacobian_ = std::make_unique(simulator_()); + diagMatAddress_.resize(numCells); + jacobian_->reserve(sparsityPattern); + for (unsigned globI = 0; globI < numCells; globI++) { + const auto& nbInfos = neighborInfo_[globI]; + diagMatAddress_[globI] = jacobian_->blockAddress(globI, globI); + for (auto& nbInfo : nbInfos) { + nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI); + } + } + + // Create full domain + fullDomain_.cells.resize(numCells); + std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0); + } + + /*! + * \brief Linearize the non-linear system of equation, computing residual and its derivatives + * + * Calculates TPSA equation terms for the AD residual, which in turn is used to get the derivates for the Jacobian + */ + template + void linearize_(const SubDomainType& domain) + { + OPM_TIMEBLOCK(linearizeTPSA); + + // Extract misc. variables used in linearization + const auto& flowModel = flowModel_(); + const auto& geoMechModel = geoMechModel_(); + auto& problem = problem_(); + const unsigned int numCells = domain.cells.size(); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + // Loop over cells in the domain and compute local residual and jacobian + for (unsigned ii = 0; ii < numCells; ++ii) { + OPM_TIMEBLOCK_LOCAL(linearizationForEachCellTPSA, Subsystem::Assembly); + + const unsigned globI = domain.cells[ii]; + const auto& nbInfos = neighborInfo_[globI]; + VectorBlock res(0.0); + MatrixBlock bMat(0.0); + ADVectorBlock adres(0.0); + const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0); + + // /// + // Face term + // /// + { + OPM_TIMEBLOCK_LOCAL(faceCalculationForEachCellTPSA, Subsystem::Assembly); + + // Loop over neighboring cells + for (const auto& nbInfo : nbInfos) { + OPM_TIMEBLOCK_LOCAL(calculationForEachFaceTPSA, Subsystem::Assembly); + + // Reset local residual and Jacobian + res = 0.0; + bMat = 0.0; + adres = 0.0; + + // Neighbor information + const unsigned globJ = nbInfo.neighbor; + assert(globJ != globI); + const MaterialState& materialStateEx = geoMechModel.materialState(globJ, /*timeIdx=*/0); + + // Compute local face term + LocalResidual::computeFaceTerm(adres, + materialStateIn, + materialStateEx, + problem, + globI, + globJ); + adres *= nbInfo.faceArea; + + // Extract residual and sub-block Jacobian entries from computed AD residual + setResAndJacobi(res, bMat, adres); + + // Insert into global residual + residual_[globI] += res; + + // Insert contribution to (globI, globI) sub-block + // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); + *diagMatAddress_[globI] += bMat; + + // Insert contribution to (globJ, globI) sub-block + // Note: since LocalResidual::computeFaceTerm have Evaluation on globI primary variables, it is + // natural to insert Jacobian entries for (globJ, globI) here, since we only need to flip the signs + // in the calculated face terms + // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat); + bMat *= -1.0; + *nbInfo.matBlockAddress += bMat; + } + } + + // /// + // Volume term + // // + adres = 0.0; + { + OPM_TIMEBLOCK_LOCAL(computeVolumeTerm, Subsystem::Assembly); + + // Compute local volume term + LocalResidual::computeVolumeTerm(adres, + materialStateIn, + problem, + globI); + } + const double volume = flowModel.dofTotalVolume(globI); + adres *= volume; + + // Extract residual and sub-block Jacobian entries from computed AD residual + setResAndJacobi(res, bMat, adres); + + // Insert in global residual + residual_[globI] += res; + + // Insert contribution to (globI, globI) sub-block + // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); + *diagMatAddress_[globI] += bMat; + + // /// + // Source term + // /// + res = 0.0; + bMat = 0.0; + adres = 0.0; + + // Compute local source term + LocalResidual::computeSourceTerm(adres, + problem, + globI, + 0); + adres *= -volume; + + // Extract residual and sub-block Jacobian entries from computed AD residual + setResAndJacobi(res, bMat, adres); + + // Insert into global residual + residual_[globI] += res; + + // Insert contribution to (globI, globI) sub-block + // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); + *diagMatAddress_[globI] += bMat; + } // globI loop + + // /// + // Boundary term + // /// + for (const auto& bdyInfo : boundaryInfo_) { + VectorBlock res(0.0); + MatrixBlock bMat(0.0); + ADVectorBlock adres(0.0); + const unsigned globI = bdyInfo.cell; + const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0); + + // Compute local boundary condition + LocalResidual::computeBoundaryTerm(adres, + materialStateIn, + bdyInfo.bcdata, + problem, + globI); + adres *= bdyInfo.bcdata.faceArea; + + // Extract residual and sub-block Jacobian entries from computed AD residual + setResAndJacobi(res, bMat, adres); + + // Insert in global residual + residual_[globI] += res; + + // Insert contribution to (globI, globI) sub-block + // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); + *diagMatAddress_[globI] += bMat; + } + } + + /*! + * \brief Create residual and Jacobian, and reset both + */ + void initFirstIteration_() + { + // initialize the BCRS matrix for the Jacobian of the residual function + createMatrix_(); + + // initialize the Jacobian matrix and the vector for the residual function + residual_.resize(flowModel_().numTotalDof()); + resetSystem_(); + } + + /*! + * \brief Reset Jacobian matrix and residuals to zero + */ + void resetSystem_() + { + // Set residual vector entries to zero + residual_ = 0.0; + + // Set all Jacobian matrix entries to zero + jacobian_->clear(); + } + + // /// + // Private get functions + // /// + /*! + * \brief Return simulator object + * + * \returns Reference to simulator + */ + Simulator& simulator_() + { return *simulatorPtr_; } + + /*! + * \brief Return simulator object + * + * \returns Reference to simulator + */ + const Simulator& simulator_() const + { return *simulatorPtr_; } + + /*! + * \brief Return problem object + * + * \returns Reference to problem + */ + Problem& problem_() + { return simulator_().problem(); } + + /*! + * \brief Return problem object + * + * \returns Reference to problem + */ + const Problem& problem_() const + { return simulator_().problem(); } + + /*! + * \brief Return Flow model object + * + * \returns Reference to Flow model + */ + FlowModel& flowModel_() + { return simulator_().model(); } + + /*! + * \brief Return Flow model object + * + * \returns Reference to Flow model + */ + const FlowModel& flowModel_() const + { return simulator_().model(); } + + /*! + * \brief Return TPSA model object + * + * \returns Reference to geomechanics model + */ + GeomechModel& geoMechModel_() + { return problem_().geoMechModel(); } + + /*! + * \brief Return TPSA model object + * + * \returns Reference to geomechanics model + */ + const GeomechModel& geoMechModel_() const + { return problem_().geoMechModel(); } + + /*! + * \brief Return grid view + * + * \returns Reference to grid view + */ + const GridView& gridView_() const + { return problem_().gridView(); } + + Simulator* simulatorPtr_{}; + LinearizationType linearizationType_{}; + + std::vector diagMatAddress_{}; + std::unique_ptr jacobian_{}; + GlobalEqVector residual_; + + // + // Helper structs + // + /*! + * \brief Neighbor information for flux terms + */ + struct NeighborInfo + { + unsigned int neighbor; + double faceArea; + MatrixBlock* matBlockAddress; + }; + SparseTable neighborInfo_{}; + + /*! + * \brief Information for boundary conditions + */ + struct BoundaryConditionData + { + BCMECHType type; + std::vector displacement; + unsigned boundaryFaceIndex; + double faceArea; + }; + + /*! + * \brief Boundary information + */ + struct BoundaryInfo + { + unsigned int cell; + int dir; + unsigned int bfIndex; + BoundaryConditionData bcdata; + }; + std::vector boundaryInfo_; + + /*! + * \brief Inforamtion of the domain to linearize + */ + struct FullDomain + { + std::vector cells; + std::vector interior; + }; + FullDomain fullDomain_; +}; // class TpsaLinearizer + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtktpsamodule.hpp b/opm/models/io/vtktpsamodule.hpp new file mode 100644 index 0000000000..1004679474 --- /dev/null +++ b/opm/models/io/vtktpsamodule.hpp @@ -0,0 +1,198 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef VTK_TPSA_MODULE_HPP +#define VTK_TPSA_MODULE_HPP + +#include + +#include +#include +#include +#include +#include +#include + + +namespace Opm { + +/*! +* \ingroup Vtk +* +* \brief VTK output module for TPSA quantities +*/ +template +class VtkTpsaModule : public BaseOutputModule +{ + using ParentType = BaseOutputModule; + + using GridView = GetPropType; + using Simulator = GetPropType; + using ElementContext = GetPropType; + + static constexpr auto vtkFormat = getPropValue(); + using VtkMultiWriter = Opm::VtkMultiWriter; + + enum { enableMech = getPropValue() }; + + using BufferType = typename ParentType::BufferType; + using ScalarBuffer = typename ParentType::ScalarBuffer; + using VectorBuffer = typename ParentType::VectorBuffer; + +public: + /*! + * \brief Constructor + * + * \param simulator Simulator object + */ + explicit VtkTpsaModule(const Simulator& simulator) + : ParentType(simulator) + { + // Read runtime parameters + if constexpr(enableMech) { + params_.read(); + } + } + + /*! + * \brief Register runtime parameters + */ + static void registerParameters() + { + if constexpr(enableMech) { + VtkTpsaParams::registerParameters(); + } + } + + /*! + * \brief Allocate memory for output + */ + void allocBuffers() override + { + // Check if vtk output is enabled + if (!Parameters::Get()) { + return; + } + + // Allocate buffers + if constexpr(enableMech) { + // Displacement + if (params_.displacementOutput_) { + this->resizeVectorBuffer_(displacement_, BufferType::Dof); + } + + // Rotation + if (params_.rotationOutput_) { + this->resizeVectorBuffer_(rotation_, BufferType::Dof); + } + + // Solid pressure + if (params_.solidPressureOutput_) { + this->resizeScalarBuffer_(solidPres_, BufferType::Dof); + } + } + } + + /*! + * \brief Assign quantities to output buffers + * + * \param elemCtx Element context object + */ + void processElement(const ElementContext& elemCtx) override + { + // Check if vtk output is enabled + if (!Parameters::Get()) { + return; + } + + if constexpr(enableMech) { + // Assign quantities from material state + const auto& problem = elemCtx.problem(); + const auto& geoMechModel = problem.geoMechModel(); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + const auto& materialState = geoMechModel.materialState(globalDofIdx, /*timeIdx=*/0); + + // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2) + for (int dirIdx = 0; dirIdx < 3; ++dirIdx) { + // Displacement + if (params_.displacementOutput_) { + displacement_[globalDofIdx][dirIdx] = scalarValue(materialState.displacement(dirIdx)); + } + + // Rotation + if (params_.rotationOutput_) { + rotation_[globalDofIdx][dirIdx] = scalarValue(materialState.rotation(dirIdx)); + } + } + + // Solid pressure + if (params_.solidPressureOutput_) { + solidPres_[globalDofIdx] = scalarValue(materialState.solidPressure()); + } + } + } + } + + /*! + * \brief Add buffers to VTK writer + * + * \param baseWriter VTK output writer object + */ + void commitBuffers(BaseOutputWriter& baseWriter) override + { + // Check if writer exists or mechanics output is enabled + if (!dynamic_cast(&baseWriter)) { + return; + } + + if constexpr(enableMech) { + // Displacement + if (params_.displacementOutput_) { + this->commitVectorBuffer_(baseWriter, "displacement", displacement_, BufferType::Dof); + } + + // Rotation + if (params_.rotationOutput_) { + this->commitVectorBuffer_(baseWriter, "rotation", rotation_, BufferType::Dof); + } + + // Solid pressure + if (params_.solidPressureOutput_) { + this->commitScalarBuffer_(baseWriter, "solid_pressure", solidPres_, BufferType::Dof); + } + } + } + +private: + VtkTpsaParams params_{}; + + VectorBuffer displacement_{}; + VectorBuffer rotation_{}; + ScalarBuffer solidPres_{}; +}; // class VtkTpsaModule + +} // namespace Opm + +#endif diff --git a/opm/models/io/vtktpsaparams.cpp b/opm/models/io/vtktpsaparams.cpp new file mode 100644 index 0000000000..de1a5a54cf --- /dev/null +++ b/opm/models/io/vtktpsaparams.cpp @@ -0,0 +1,54 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include +#include + + +namespace Opm { + +/*! +* \brief Register runtime parameters +*/ +void VtkTpsaParams::registerParameters() +{ + Parameters::Register + ("Include displacement in VTK output files"); + Parameters::Register + ("Include rotation in VTK output files"); + Parameters::Register + ("Include solid pressure in VTK output files"); +} + +/*! +* \brief Read runtime parameters +*/ +void VtkTpsaParams::read() +{ + displacementOutput_ = Parameters::Get(); + rotationOutput_ = Parameters::Get(); + solidPressureOutput_ = Parameters::Get(); +} + +} // namespace Opm diff --git a/opm/models/io/vtktpsaparams.hpp b/opm/models/io/vtktpsaparams.hpp new file mode 100644 index 0000000000..344dadfd86 --- /dev/null +++ b/opm/models/io/vtktpsaparams.hpp @@ -0,0 +1,55 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef VTK_TPSA_PARAMS_HPP +#define VTK_TPSA_PARAMS_HPP + + +namespace Opm::Parameters { + +// Default output quantities +struct VtkWriteDisplacement { static constexpr bool value = true; }; +struct VtkWriteRotation { static constexpr bool value = true; }; +struct VtkWriteSolidPressure { static constexpr bool value = true; }; + +} // namespace Opm::Parameters + +namespace Opm { + +/*! +* \brief Parameters for VtkTpsaOutputModule +*/ +struct VtkTpsaParams +{ + static void registerParameters(); + void read(); + + bool displacementOutput_; + bool rotationOutput_; + bool solidPressureOutput_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/tpsa/elasticityindices.hpp b/opm/models/tpsa/elasticityindices.hpp new file mode 100644 index 0000000000..8ced845c51 --- /dev/null +++ b/opm/models/tpsa/elasticityindices.hpp @@ -0,0 +1,64 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef ELASTICITY_INDICES_HPP +#define ELASTICITY_INDICES_HPP + + +namespace Opm { + +template +struct ElasticityIndices +{ + // /// + // Primary variable indices + // /// + // Starting index of displacement vector + static constexpr int disp0Idx = PVOffset + 0; + + // Starting index of rotation vector + static constexpr int rot0Idx = disp0Idx + 3; + + // Index of solid pressure + static constexpr int solidPres0Idx = rot0Idx + 3; + + // /// + // Equation indices + // /// + // Starting index of all equations, in particular the stress equations + static constexpr int conti0EqIdx = PVOffset + 0; + + // Starting index of rotation equations + static constexpr int contiRotEqIdx = conti0EqIdx + 3; + + // Starting index of solid pressure equation + static constexpr int contiSolidPresEqIdx = contiRotEqIdx + 3; + + // Total number of equations/primary variables + static constexpr int numEq = 7; +}; // struct ElasticityIndices + +} // namespace Opm + +#endif diff --git a/opm/models/tpsa/elasticitylocalresidualtpsa.hpp b/opm/models/tpsa/elasticitylocalresidualtpsa.hpp new file mode 100644 index 0000000000..bb408a0fdc --- /dev/null +++ b/opm/models/tpsa/elasticitylocalresidualtpsa.hpp @@ -0,0 +1,426 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef ELASTICITY_LOCAL_RESIDUAL_TPSA_HPP +#define ELASTICITY_LOCAL_RESIDUAL_TPSA_HPP + +#include + +#include + +#include +#include + +#include + +#include + +#include + + +namespace Opm { + +/*! +* \brief Calculation of (linear) elasticity model terms for the residual +* +* The linearized Biot model is solved where it is assumed that solid mechanics are governed by Hooke's law and +* conservation of linear momentum: +* +* \f{eqnarray*}{ +* \sigma &=& 2\mu\epsilon(u) + \lambda(\nabla\cdot u)I, +* -\nabla\cdot(\sigma - \alpha(p_f - p_0)I) &=& f_u, +* \f} +* +* where \f&\sigma\f& is Cauchy stress tensor, \f&u\f& is displacement, \f&\epsilon(\cdot)\f& is thesymmetric gradient, +* \f&\mu\f& and \f&\lambda\f& are Lame's first (aka shear modulus) and second parameters, \f&\alpha\f& is the +* Biot-Willis parameter, \f&(p_f-p_0)\f& is fluid pressure difference wrt hydrostatic, and \f&f_u\f& are body forces. +* +* The equations are discretized using two-point stress approximation following Boon et al. (2025), Solving Biot +* poroelasticity by coupling OPM Flow with the two-point stress approximation finite volume method, arXiv:2510.23432v1. +* +* The resulting equations contain a volume term where only single-cell variables are used; face terms where variables +* across cell faces are calculated; boundary terms, similar to face terms, but cell faces are at the boundary; and +* source terms where coupling and potential body forces are calculated. +*/ +template +class ElasticityLocalResidual +{ + using Indices = GetPropType; + using Scalar = GetPropType; + using Evaluation = GetPropType; + using Problem = GetPropType; + + using MaterialState = MaterialStateTPSA; + using Toolbox = MathToolbox; + + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { contiRotEqIdx = Indices::contiRotEqIdx }; + enum { contiSolidPresEqIdx = Indices::contiSolidPresEqIdx }; + enum { numEq = getPropValue() }; + +public: + /*! + * \brief Calculate volume terms in TPSA formulation + * + * \param volTerm Volume term vector + * \param materialState Material state container + * \param problem Flow problem + * \param globalIndex Cell index + * + * Material state, problem input and global index here might/should be merged in an "IntensiveQuantity" container as + * in BlackOilLocalResidualTPFA + */ + template + static void computeVolumeTerm(Dune::FieldVector& volTerm, + const MaterialState& materialState, + const Problem& problem, + const unsigned globalIndex) + { + // Reset volume terms + volTerm = 0.0; + + // Rotation equations (one per direction) + const Scalar sModulus = problem.shearModulus(globalIndex); + for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) { + volTerm[contiRotEqIdx + dirIdx] += + Toolbox::template decay(materialState.rotation(dirIdx)) + / sModulus; + } + + // Solid pressure equation + const Scalar lame = problem.lame(globalIndex); + volTerm[contiSolidPresEqIdx] += + Toolbox::template decay(materialState.solidPressure()) + / lame; + } + + /*! + * \brief Calculate terms across cell faces in TPSA formulation + * + * \param faceTerm Face term vector + * \param materialStateIn Material state container of inside cell + * \param materialStateEx Material state container of outside cel + * \param problem Flow problem + * \param globalIndexIn Inside cell index + * \param globalIndexEx Outside cell index + * + * Material state, problem input and global index here might/should be merged in "IntensiveQuantity" and + * "NeighborInfo" containers as in BlackOilLocalResidualTPFA + */ + static void computeFaceTerm(Dune::FieldVector& faceTerm, + const MaterialState& materialStateIn, + const MaterialState& materialStateEx, + Problem& problem, + const unsigned globalIndexIn, + const unsigned globalIndexEx) + { + // Reset face terms + faceTerm = 0.0; + + // Extract some face properties + const Scalar weightAvgIn = problem.weightAverage(globalIndexIn, globalIndexEx); + const Scalar weightAvgEx = problem.weightAverage(globalIndexEx, globalIndexIn); + const Scalar weightProd = problem.weightProduct(globalIndexIn, globalIndexEx); + const Scalar normDist = problem.normalDistance(globalIndexIn, globalIndexEx); + const auto& faceNormal = problem.cellFaceNormal(globalIndexIn, globalIndexEx); + + // Effective shear modulus + const Scalar sModulusIn = problem.shearModulus(globalIndexIn); + const Scalar sModulusEx = problem.shearModulus(globalIndexEx); + const Scalar eff_sModulus = weightAvgIn * sModulusIn + weightAvgEx * sModulusEx; + + // Distance ratio + const Scalar distRatio = 0.5 * weightProd / normDist; + + // Solid pressures + const Evaluation& solidPIn = materialStateIn.solidPressure(); + const Scalar solidPEx = decay(materialStateEx.solidPressure()); + + // /// + // Solid pressure equation (direction-independent equation) + // /// + faceTerm[contiSolidPresEqIdx] += + distRatio * eff_sModulus * (solidPIn - solidPEx); + + // /// + // Displacement, rotation and solid pressure (directional-dependent) equations + // /// + // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product. + // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1 + auto modNeg = [](int i) { return ((i % 3) + 3) % 3; }; + + // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2) + for (int dirIdx = 0; dirIdx < 3; ++dirIdx) { + // Direction indices in cross-product + unsigned dirIdxNeg = modNeg(dirIdx - 1); + unsigned dirIdxPos = modNeg(dirIdx + 1); + + // Displacement equation + const Scalar faceNormalDir = faceNormal[dirIdx]; + const Scalar faceNormalNeg = faceNormal[dirIdxNeg]; + const Scalar faceNormalPos = faceNormal[dirIdxPos]; + + const Evaluation& dispIn = materialStateIn.displacement(dirIdx); + const Scalar dispEx = decay(materialStateEx.displacement(dirIdx)); + + const Evaluation& rotInNeg = materialStateIn.rotation(dirIdxNeg); + const Evaluation& rotInPos = materialStateIn.rotation(dirIdxPos); + const Scalar rotExNeg = decay(materialStateEx.rotation(dirIdxNeg)); + const Scalar rotExPos = decay(materialStateEx.rotation(dirIdxPos)); + + faceTerm[conti0EqIdx + dirIdx] += + 2.0 * (eff_sModulus / normDist) * (dispIn - dispEx) + - weightAvgIn * (faceNormalNeg * rotInPos - faceNormalPos * rotInNeg) + - weightAvgEx * (faceNormalNeg * rotExPos - faceNormalPos * rotExNeg) + - faceNormalDir * (weightAvgIn * solidPIn + weightAvgEx * solidPEx); + + // Rotation equation + const Evaluation& dispInNeg = materialStateIn.displacement(dirIdxNeg); + const Evaluation& dispInPos = materialStateIn.displacement(dirIdxPos); + const Scalar dispExNeg = decay(materialStateEx.displacement(dirIdxNeg)); + const Scalar dispExPos = decay(materialStateEx.displacement(dirIdxPos)); + + faceTerm[contiRotEqIdx + dirIdx] += + - weightAvgEx * (faceNormalNeg * dispInPos - faceNormalPos * dispInNeg) + - weightAvgIn * (faceNormalNeg * dispExPos - faceNormalPos * dispExNeg); + + // Solid pressure equation + faceTerm[contiSolidPresEqIdx] += + - faceNormalDir * (weightAvgEx * dispIn + weightAvgIn * dispEx); + } + } + + /*! + * \brief Calculate boundary conditions in TPSA formulation given by BCCON/BCPROP + * + * \param bndryTerm Boundary term vector + * \param materialState Material state container + * \param bdyInfo Boundary condition info container + * \param problem Flow problem + * \param globalIndex Cell index + */ + template + static void computeBoundaryTerm(Dune::FieldVector& bndryTerm, + const MaterialState& materialState, + const BoundaryConditionData& bdyInfo, + Problem& problem, + unsigned globalIndex) + { + // Switch between possible boundary conditions + switch (bdyInfo.type) { + // OBS: NONE is interpreted as FIXED with zero displacement + case BCMECHType::NONE: + computeBoundaryTermFixed(bndryTerm, + materialState, + bdyInfo, + problem, + globalIndex); + break; + case BCMECHType::FIXED: + throw std::runtime_error("BCTYPE FIXED has not been implemented in TPSA"); + case BCMECHType::FREE: + computeBoundaryTermFree(bndryTerm, + materialState, + bdyInfo, + problem, + globalIndex); + break; + default: + throw std::logic_error("Unknown boundary condition type " + + std::to_string(static_cast(bdyInfo.type)) + + " in computeBoundaryFlux()." ); + } + } + + /*! + * \brief Calculate fixed displacement boundary condition in TPSA formulation + * + * \param bndryTerm Boundary term vector + * \param materialState Material state container + * \param bdyInfo Boundary condition info container + * \param problem Flow problem + * \param globalIndex Cell index + * + * \note BCMECHType::NONE is a implemented as a specialization of BCMECHType::FIXED with displacement equal + * to zero on the boundary face + */ + template + static void computeBoundaryTermFixed(Dune::FieldVector& bndryTerm, + const MaterialState& materialState, + const BoundaryConditionData& bdyInfo, + Problem& problem, + unsigned globalIndex) + { + // !!!! + // Only BCMECHType::NONE, where we have zero displacement on boundary face, have been implemented! + // !!!! + + // Reset bondary term + bndryTerm = 0.0; + + // Extract some face properties + const unsigned bfIdx = bdyInfo.boundaryFaceIndex; + const Scalar weightAvg = problem.weightAverageBoundary(globalIndex, bfIdx); + const Scalar normDist = problem.normalDistanceBoundary(globalIndex, bfIdx); + const auto& faceNormal = problem.cellFaceNormalBoundary(globalIndex, bfIdx); + + // Effective shear modulus (= cell shear modulus) + const Scalar eff_sModulus = problem.shearModulus(globalIndex); + + // Solid pressure + const Evaluation& solidP = materialState.solidPressure(); + + // /// + // Displacement equation + // /// + // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product. + // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1 + auto modNeg = [](int i) { return ((i % 3) + 3) % 3; }; + + // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2) + for (int dirIdx = 0; dirIdx < 3; ++dirIdx) { + // Direction indices in cross-product + unsigned dirIdxNeg = modNeg(dirIdx - 1); + unsigned dirIdxPos = modNeg(dirIdx + 1); + + // Displacement equation + const Scalar faceNormalDir = faceNormal[dirIdx]; + const Scalar faceNormalNeg = faceNormal[dirIdxNeg]; + const Scalar faceNormalPos = faceNormal[dirIdxPos]; + + const Evaluation& disp = materialState.displacement(dirIdx); + + const Evaluation& rotNeg = materialState.rotation(dirIdxNeg); + const Evaluation& rotPos = materialState.rotation(dirIdxPos); + + bndryTerm[conti0EqIdx + dirIdx] += + 2.0 * (eff_sModulus / normDist) * disp + - weightAvg * (faceNormalNeg * rotPos - faceNormalPos * rotNeg) + - faceNormalDir * weightAvg * solidP; + } + } + + /*! + * \brief Calculate free (or zero traction) boundary condition in TPSA formulation + * + * \param bndryTerm Boundary term vector + * \param materialState Material state container + * \param bdyInfo Boundary condition info container + * \param problem Flow problem + * \param globalIndex Cell index + * + * \note Free, or zero traction, BC is equivalent of having a spring at infinity where we have assumed all (primary) + * variables and parameters (e.g., shear modulus) are zero. + */ + template + static void computeBoundaryTermFree(Dune::FieldVector& bndryTerm, + const MaterialState& materialState, + const BoundaryConditionData& bdyInfo, + Problem& problem, + unsigned globalIndex) + { + // Reset bondary term + bndryTerm = 0.0; + + // Face properties + const unsigned bfIdx = bdyInfo.boundaryFaceIndex; + const Scalar weightAvg = 1.0; + const Scalar normDist = problem.normalDistanceBoundary(globalIndex, bfIdx); + const auto& faceNormal = problem.cellFaceNormalBoundary(globalIndex, bfIdx); + + const Scalar sModulus = problem.shearModulus(globalIndex); + + // /// + // Solid pressure equation (direction-independent equation) + // /// + const Evaluation& solidP = materialState.solidPressure(); + bndryTerm[contiSolidPresEqIdx] += + 0.5 * (normDist / sModulus) * solidP; + + // /// + // Rotation and solid pressure (directional-dependent) equations + // /// + // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product. + // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1 + auto modNeg = [](int i) { return ((i % 3) + 3) % 3; }; + + // Pre-compute dot product for rotation equation + Evaluation dotProd = 0; + for (int dirIdx = 0; dirIdx < 3; ++dirIdx) { + dotProd += faceNormal[dirIdx] * materialState.rotation(dirIdx); + } + + // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2) + for (int dirIdx = 0; dirIdx < 3; ++dirIdx) { + // Direction indices in cross-product + unsigned dirIdxNeg = modNeg(dirIdx - 1); + unsigned dirIdxPos = modNeg(dirIdx + 1); + + // Rotation equation + const Scalar faceNormalDir = faceNormal[dirIdx]; + const Scalar faceNormalNeg = faceNormal[dirIdxNeg]; + const Scalar faceNormalPos = faceNormal[dirIdxPos]; + + const Evaluation& dispNeg = materialState.displacement(dirIdxNeg); + const Evaluation& dispPos = materialState.displacement(dirIdxPos); + + const Evaluation& rot = materialState.rotation(dirIdx); + + bndryTerm[contiRotEqIdx + dirIdx] += + - weightAvg * (faceNormalNeg * dispPos - faceNormalPos * dispNeg) + + 0.5 * (normDist / sModulus) * (dotProd * faceNormalDir - rot); + + // Solid pressure (directional-dependent) equation + const Evaluation& disp = materialState.displacement(dirIdx); + + bndryTerm[contiSolidPresEqIdx] += + - faceNormalDir * weightAvg * disp; + } + } + + /*! + * \brief Calculate source term in TPSA formulation + * + * \param sourceTerm Source term vector + * \param problem Flow problem + * \param globalSpaceIdex Cell index + * \param timeIdx Time index + */ + static void computeSourceTerm(Dune::FieldVector& sourceTerm, + Problem& problem, + unsigned globalSpaceIdex, + unsigned timeIdx) + { + // Reset source terms + sourceTerm = 0.0; + + // Get source term from problem + // NOTE: separate source function than Flow source(...)! + problem.tpsaSource(sourceTerm, globalSpaceIdex, timeIdx); + } +}; // class ElasticityLocalResidual + +} // namespace Opm + +#endif diff --git a/opm/models/tpsa/elasticityprimaryvariables.hpp b/opm/models/tpsa/elasticityprimaryvariables.hpp new file mode 100644 index 0000000000..944b45d322 --- /dev/null +++ b/opm/models/tpsa/elasticityprimaryvariables.hpp @@ -0,0 +1,153 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef ELASTICITY_PRIMARY_VARIABLES_HPP +#define ELASTICITY_PRIMARY_VARIABLES_HPP + +#include + +#include +#include + +#include +#include +#include + +#include + + +namespace Opm { + +/*! +* \brief Primary variables in (linear) elasticity equations +* +* Primary variables are: +* \li Displacement (3D vector) +* \li Rotation (3D vector) -> variable to express rotations in stress tensor +* \li Solid pressure (single scalar) -> volumetric change from mechanics +*/ +template +class ElasticityPrimaryVariables + : public Dune::FieldVector, + getPropValue()> +{ + using Evaluation = GetPropType; + using Indices = GetPropType; + using Scalar = GetPropType; + + using Toolbox = Opm::MathToolbox; + + enum { disp0Idx = Indices::disp0Idx }; + enum { rot0Idx = Indices::rot0Idx }; + enum { solidPres0Idx = Indices::solidPres0Idx }; + + enum { numEq = getPropValue() }; + + using ParentType = Dune::FieldVector; + +public: + /*! + * \brief Constructor + */ + ElasticityPrimaryVariables() : ParentType() + { + Valgrind::SetUndefined(*this); + } + + /*! + * \brief Default copy constructor + */ + ElasticityPrimaryVariables(const ElasticityPrimaryVariables& value) = default; + + /*! + * \brief Default assignment constructor + */ + ElasticityPrimaryVariables& operator=(const ElasticityPrimaryVariables& value) = default; + + using ParentType::operator=; //!< Import base class assignment operators. + + /*! + * \brief Return primary variable in Evaluation type + * + * \param varIdx Primary variable index + * \param timeIdx Time index + * \param linearizationType Type of linearization + * \returns Primary variable as Evalutation type + * + * \note + * \li Automatic differentiation: returns value + derivative + * \li Finite differences: returns value only + */ + Evaluation makeEvaluation(unsigned varIdx, + unsigned timeIdx, + Opm::LinearizationType linearizationType = LinearizationType()) const + { + // Finite difference + if constexpr (std::is_same_v) { + return (*this)[varIdx]; + } + // Automatic differentiation + else { + if (timeIdx == linearizationType.time) { + return Toolbox::createVariable((*this)[varIdx], varIdx); + } + else { + return Toolbox::createConstant((*this)[varIdx]); + } + } + } + + /*! + * \brief Assign primary variables from a material state container + * + * \param materialState Material state container + */ + template + void assignNaive(const MaterialState& materialState) + { + // Reset primary variables vector + (*this) = 0.0; + + // Assign displacement and rotation vectors + for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) { + (*this)[disp0Idx + dirIdx] = Opm::getValue(materialState.displacement(dirIdx)); + (*this)[rot0Idx + dirIdx] = Opm::getValue(materialState.rotation(dirIdx)); + } + + // Assign solid pressure + (*this)[solidPres0Idx] = Opm::getValue(materialState.solidPressure()); + } + + /*! + * \brief Instruct Valgrind to check the definedness of all attributes of this class + */ + void checkDefined() const + { + Valgrind::CheckDefined(*static_cast(this)); + } +}; // class ElasticityPrimaryVariables + +} // namespace Opm + +#endif diff --git a/opm/models/tpsa/tpsabaseproperties.hpp b/opm/models/tpsa/tpsabaseproperties.hpp new file mode 100644 index 0000000000..12abd3e3fb --- /dev/null +++ b/opm/models/tpsa/tpsabaseproperties.hpp @@ -0,0 +1,77 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_BASE_PROPERTIES_HPP +#define TPSA_BASE_PROPERTIES_HPP + +#include + + +namespace Opm::Properties { + +template +struct IndicesTPSA { using type = UndefinedProperty; }; + +template +struct NumEqTPSA { using type = UndefinedProperty; }; + +template +struct LinearizerTPSA { using type = UndefinedProperty; }; + +template +struct EvaluationTPSA { using type = UndefinedProperty; }; + +template +struct EqVectorTPSA { using type = UndefinedProperty; }; + +template +struct GlobalEqVectorTPSA { using type = UndefinedProperty; }; + +template +struct PrimaryVariablesTPSA { using type = UndefinedProperty; }; + +template +struct SolutionVectorTPSA { using type = UndefinedProperty; }; + +template +struct SolutionHistorySizeTPSA { using type = UndefinedProperty; }; + +template +struct ModelTPSA { using type = UndefinedProperty; }; + +template +struct LocalResidualTPSA { using type = UndefinedProperty; }; + +template +struct SparseMatrixAdapterTPSA { using type = UndefinedProperty; }; + +template +struct NewtonMethodTPSA { using type = UndefinedProperty; }; + +template +struct LinearSolverBackendTPSA { using type = UndefinedProperty; }; + +} // namespace Opm::Properties + +#endif diff --git a/opm/models/tpsa/tpsamodel.hpp b/opm/models/tpsa/tpsamodel.hpp new file mode 100644 index 0000000000..c2272927cc --- /dev/null +++ b/opm/models/tpsa/tpsamodel.hpp @@ -0,0 +1,582 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_MODEL_HPP +#define TPSA_MODEL_HPP + +#include + +#include + +#include + +#include +#include +#include + +#include +#include + + +namespace Opm { + +/*! +* \brief TPSA geomechanics model +* +* Solves the (linear) elasticity equations using Newton method. +*/ +template +class TpsaModel +{ + using DofMapper = GetPropType; + using Evaluation = GetPropType; + using GlobalEqVector = GetPropType; + using GridView = GetPropType; + using Indices = GetPropType; + using Linearizer = GetPropType; + using NewtonMethod = GetPropType; + using PrimaryVariables = GetPropType; + using Simulator = GetPropType; + using SolutionVector = GetPropType; + using Scalar = GetPropType; + + enum { dimWorld = GridView::dimensionworld }; + enum { historySize = getPropValue() }; + enum { numEq = getPropValue() }; + + enum { disp0Idx = Indices::disp0Idx }; + enum { rot0Idx = Indices::rot0Idx }; + enum { solidPres0Idx = Indices::solidPres0Idx }; + + using MaterialState = MaterialStateTPSA; + + using DimVector = Dune::FieldVector; + using SymTensor = Dune::FieldVector; + +public: + /*! + * \brief Small block vector wrapper class for model solutions + */ + class TpsaBlockVectorWrapper + { + protected: + SolutionVector blockVector_; + public: + /*! + * \brief Constructor + * + * \param + * \param size Size of block vector + */ + TpsaBlockVectorWrapper(const std::string&, const std::size_t size) + : blockVector_(size) + {} + + /*! + * \brief Default constructor + */ + TpsaBlockVectorWrapper() = default; + + /*! + * \brief Test function for serialization + */ + static TpsaBlockVectorWrapper serializationTestObject() + { + TpsaBlockVectorWrapper result("dummy", 3); + result.blockVector_[0] = 1.0; + result.blockVector_[1] = 2.0; + result.blockVector_[2] = 3.0; + + return result; + } + + /*! + * \brief Get reference of block vector + * + * \returns Block vector + */ + SolutionVector& blockVector() + { return blockVector_; } + + /*! + * \brief Get const reference of block vector + * + * \returns Block vector + */ + const SolutionVector& blockVector() const + { return blockVector_; } + + /*! + * \brief Check if incoming block vector is the same as current + * + * \param wrapper Block vector to check + * \returns Boolean indicating if wrapper is equal to current block vector + */ + bool operator==(const TpsaBlockVectorWrapper& wrapper) const + { + return std::equal(this->blockVector_.begin(), this->blockVector_.end(), + wrapper.blockVector_.begin(), wrapper.blockVector_.end()); + } + + /*! + * \brief Serializing operation + * + * \param serializer Reference to serializer object + */ + template + void serializeOp(Serializer& serializer) + { + serializer(blockVector_); + } + }; + + // /// + // Public functions + // /// + /*! + * \brief Constructor + * + * \param simulator Simulator object + */ + explicit TpsaModel(Simulator& simulator) + : linearizer_(std::make_unique()) + , newtonMethod_(simulator) + , simulator_(simulator) + , element_chunks_(simulator.gridView(), Dune::Partitions::all, ThreadManager::maxThreads()) + { + // Initialize equation weights to 1.0 + eqWeights_.resize(numEq, 1.0); + + // Initialize historic solution vectors + // OBS: need at least history size = 2, due to time-derivative of solid-pressure in Flow coupling term + const std::size_t numDof = simulator_.model().numGridDof(); + for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) { + solution_[timeIdx] = std::make_unique("solution", numDof); + } + } + + /*! + * \brief Initialize TPSA model + */ + void finishInit() + { + // Initialize the linearizer + linearizer_->init(simulator_); + + // Resize material state vector + resizeMaterialState_(); + } + + /*! + * \brief Register runtime parameters + */ + static void registerParameters() + { + // Newton method parameters + NewtonMethod::registerParameters(); + } + + /*! + * \brief Prepare TPSA model for coupled Flow-TPSA scheme + */ + void prepareTPSA() + { + // Update historic solution + solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0); + } + + /*! + * \brief Sync primary variables in overlapping cells + * + * Copied code from EcfvDiscretization::syncOverlap() to sync TPSA primary variables + */ + void syncOverlap() + { + // Syncronize the solution on the ghost and overlap elements + using GhostSyncHandle = GridCommHandleGhostSync; + + auto ghostSync = GhostSyncHandle(solution(/*timeIdx=*/0), + simulator_.model().dofMapper()); + + simulator_.gridView().communicate(ghostSync, + Dune::InteriorBorder_All_Interface, + Dune::ForwardCommunication); + } + + /*! + * \brief Update material state for all cells + * + * \param timeIdx Time index + * + * \note Cached material state not implemented yet! + */ + void updateMaterialState(const unsigned /*timeIdx*/) + { + // Loop over all elements chuncks and update material state from current solution + const auto& elementMapper = simulator_.model().elementMapper(); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (const auto& chunk : element_chunks_) { + for (const auto& elem : chunk) { + const unsigned globalIdx = elementMapper.index(elem); + auto& currSol = solution(/*timeIdx=*/0)[globalIdx]; + setMaterialState_(globalIdx, /*timeIdx=*/0, currSol); + } + } + } + + // /// + // Public get and set functions + // /// + /*! + * \brief Return the linearizer + * + * \returns Reference to linearizer + */ + const Linearizer& linearizer() const + { + return *linearizer_; + } + + /*! + * \brief Return the linearizer + * + * \returns Reference to linearizer + */ + Linearizer& linearizer() + { + return *linearizer_; + } + + /*! + * \brief Return the Newton method + * + * \returns Reference to Newton method object + */ + const NewtonMethod& newtonMethod() const + { + return newtonMethod_; + } + + /*! + * \brief Return the Newton method + * + * \returns Reference to Newton method object + */ + NewtonMethod& newtonMethod() + { + return newtonMethod_; + } + + /*! + * \brief Get reference to history solution vector + * + * \param timeIdx Time index + * \returns Reference to solution vector at time step timeIdx + */ + const SolutionVector& solution(unsigned timeIdx) const + { + return solution_[timeIdx]->blockVector(); + } + + /*! + * \brief Get reference to history solution vector + * + * \param timeIdx Time index + * \returns Reference to solution vector at time step timeIdx + */ + SolutionVector& solution(unsigned timeIdx) + { + return solution_[timeIdx]->blockVector(); + } + + /*! + * \brief Return number of degrees of freedom in the grid from the Flow model + * \returns Number of grid DOFs + */ + std::size_t numGridDof() const + { + return simulator_.model().numGridDof(); + } + + /*! + * \brief Return the total number of degrees of freedom + * \returns Number of grid DOFs + auxillary DOFs + */ + std::size_t numTotalDof() const + { + return numGridDof() + numAuxiliaryDof(); + } + + /*! + * \brief Return the total grid volume from the Flow model + * + * \param globalIdx Cell index + * \returns Grid volume of grid cell + */ + Scalar dofTotalVolume(unsigned globalIdx) const + { + return simulator_.model().dofTotalVolume(globalIdx); + } + + /*! + * \brief Return equation weights + * + * \param dofIdx Degree-of-freedom index + * \param eqIdx Equation index + * \returns Weight for equation + */ + Scalar eqWeight(unsigned /*dofIdx*/, unsigned eqIdx) const + { + return eqWeights_[eqIdx]; + } + + /*! + * \brief Set weights for equation + * + * \param eqIdx Equation index + * \param value Weight value to set + */ + void setEqWeight(unsigned eqIdx, Scalar value) + { + eqWeights_[eqIdx] = value; + } + + /*! + * \brief Return number of auxillary modules + * \returns No. of auxillary modules + */ + std::size_t numAuxiliaryModules() const + { + return 0; + } + + /*! + * \brief Return number of auxillary degrees of freedom + * \returns No. of auxillary DOFs + */ + std::size_t numAuxiliaryDof() const + { + return 0; + } + + /*! + * \brief Return current material state + * + * \param globalIdx Cell index + * \param timeIdx Time index + * \returns Material state container for grid cell + * + * \note Cached material state not implemented yet! + */ + const MaterialState& materialState(const unsigned globalIdx, unsigned /*timeIdx*/) const + { + return materialState_[globalIdx]; + } + + /*! + * \brief Output displacement vector + * + * \param globalIdx Cell index + * \param with_fracture Boolean to activate fracture output + * \returns Displacement vector at grid cell + * + * \note Needed in OutputBlackOilModule! + */ + DimVector disp(const unsigned globalIdx, const bool /*with_fracture*/) const + { + DimVector d; + for (std::size_t i = 0; i < 3; ++i) { + d[i] = decay(materialState_[globalIdx].displacement(i)); + } + return d; + } + + /*! + * \brief Output (del?)stress tensor + * + * \param globalIdx Cell index + * \returns (Del?)stress tensor (Voigt notation) at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + SymTensor delstress(const unsigned /*globalIdx*/) const + { + SymTensor val; + return val; + } + + /*! + * \brief Output fracture stress tensor + * + * \param globalIdx Cell index + * \returns Fracture stress tensor (Voigt notation) at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + SymTensor fractureStress(const unsigned /*globalIdx*/) const + { + SymTensor val; + return val; + } + + /*! + * \brief Output linear stress tensor + * + * \param globalIdx Cell index + * \returns Linear stress tensor (Voigt notation) at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + SymTensor linstress(const unsigned /*globalIdx*/) const + { + SymTensor val; + return val; + } + + /*! + * \brief Output stress tensor + * + * \param globalIdx Cell index + * \param with_fracture Boolean to activate fracture output + * \returns Stress tensor (Voigt notation) at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + SymTensor stress(const unsigned /*globalIdx*/, const bool /*with_fracture*/) const + { + SymTensor val; + return val; + } + + /*! + * \brief Output strain tensor + * + * \param globalIdx Cell index + * \param with_fracture Boolean to activate fracture output + * \returns Strain tensor (Voigt notation) at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + SymTensor strain(const unsigned /*globalIdx*/, const bool /*with_fracture*/) const + { + SymTensor val; + return val; + } + + /*! + * \brief Output potential forces + * + * \param globalIdx Cell index + * \returns Potential forces at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + Scalar mechPotentialForce(unsigned /*globalIdx*/) const + { + return Scalar(0.0); + } + + /*! + * \brief Output potential pressure forces + * + * \param globalIdx Cell index + * \returns Potential pressure forces at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + Scalar mechPotentialPressForce(unsigned /*globalIdx*/) const + { + return Scalar(0.0); + } + + /*! + * \brief Output potential temparature forces + * + * \param globalIdx Cell index + * \returns Potential temperature forces at grid cell + * + * \note Needed in OutputBlackOilModule, but zero for now! + */ + Scalar mechPotentialTempForce(unsigned /*globalIdx*/) const + { + return Scalar(0.0); + } + +protected: + // /// + // Protected functions + // /// + /*! + * \brief Resize material state vector + * + * \note Cached material state not implemented yet! + */ + void resizeMaterialState_() + { + const std::size_t numDof = simulator_.model().numGridDof(); + materialState_.resize(numDof); + } + +private: + // /// + // Private functions + // /// + /*! + * \brief Set material from another + * + * \param globalIdx Cell index + * \param timeIdx Time index + * \param values Primary variable to insert into material state + * + * \note Cached material state not implemented yet! + */ + void setMaterialState_(const unsigned globalIdx, const unsigned /*timeIdx*/, PrimaryVariables& values) + { + auto& dofMaterialState = materialState_[globalIdx]; + for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) { + dofMaterialState.setDisplacement(dirIdx, values.makeEvaluation(disp0Idx + dirIdx, 0)); + dofMaterialState.setRotation(dirIdx, values.makeEvaluation(rot0Idx + dirIdx, 0)); + } + dofMaterialState.setSolidPressure(values.makeEvaluation(solidPres0Idx, 0)); + } + + std::unique_ptr linearizer_; + NewtonMethod newtonMethod_; + Simulator& simulator_; + ElementChunks element_chunks_; + + std::array, historySize> solution_; + std::vector eqWeights_; + std::vector materialState_; +}; // class TpsaModel + +} // namespace Opm + + +#endif diff --git a/opm/models/tpsa/tpsanewtonmethod.hpp b/opm/models/tpsa/tpsanewtonmethod.hpp new file mode 100644 index 0000000000..387838a5c2 --- /dev/null +++ b/opm/models/tpsa/tpsanewtonmethod.hpp @@ -0,0 +1,578 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_NEWTON_METHOD_HPP +#define TPSA_NEWTON_METHOD_HPP + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace Opm { + +/*! +* \brief Newton method solving for generic TPSA model. +* +* Generates the Jacobian matrix, J(u^n) and residual vector, R(u^n), with a solution vector, u^n, at iteration n. +* Subsequently the linear system J(u^n)\delta u^n = -R(u^n) is solved to get u^{n+1} = u^n + \Delta u^n. +*/ +template +class TpsaNewtonMethod +{ + using Scalar = GetPropType; + using Simulator = GetPropType; + using Problem = GetPropType; + using Model = GetPropType; + + using SolutionVector = GetPropType; + using GlobalEqVector = GetPropType; + using PrimaryVariables = GetPropType; + using EqVector = GetPropType; + using Linearizer = GetPropType; + using LinearSolverBackend = GetPropType; + using SparseMatrixAdapter = GetPropType; + + using IstlMatrix = typename SparseMatrixAdapter::IstlMatrix; + +public: + /*! + * \brief Constructor + * + * \param simulator Simulator object + */ + explicit TpsaNewtonMethod(Simulator& simulator) + : simulator_(simulator) + , linearSolver_(simulator) + , error_(1e100) + , lastError_(1e100) + , numIterations_(0) + , numLinearizations_(0) + { + // Read runtime/default Newton parameters + params_.read(); + } + + /*! + * \brief Register all run-time parameters for the Newton method. + */ + static void registerParameters() + { + LinearSolverBackend::registerParameters(); + TpsaNewtonMethodParams::registerParameters(); + } + + /*! + * \brief Run the Newton method. + * + * \returns Bool indicating if Newton converged + */ + bool apply() + { + // Make sure all timers are prestine + prePostProcessTimer_.halt(); + linearizeTimer_.halt(); + solveTimer_.halt(); + updateTimer_.halt(); + + // Get vectors and linearizer + SolutionVector& nextSolution = model().solution(/*historyIdx=*/0); + SolutionVector currentSolution(nextSolution); + GlobalEqVector solutionUpdate(nextSolution.size()); + + Linearizer& linearizer = model().linearizer(); + + TimerGuard prePostProcessTimerGuard(prePostProcessTimer_); + + // Tell the implementation that we begin solving + prePostProcessTimer_.start(); + begin_(); + prePostProcessTimer_.stop(); + + try { + TimerGuard innerPrePostProcessTimerGuard(prePostProcessTimer_); + TimerGuard linearizeTimerGuard(linearizeTimer_); + TimerGuard updateTimerGuard(updateTimer_); + TimerGuard solveTimerGuard(solveTimer_); + + // Execute the method as long as the implementation thinks that we should do another iteration + while (proceed_()) { + // Notify the implementation that we're about to start a new iteration + prePostProcessTimer_.start(); + beginIteration_(); + prePostProcessTimer_.stop(); + + // Make the current solution to the old one + currentSolution = nextSolution; + + // Do the actual linearization + linearizeTimer_.start(); + linearizeDomain_(); + linearizeTimer_.stop(); + + // Get residual and Jacobian for convergence check and preparation of linear solver + solveTimer_.start(); + auto& residual = linearizer.residual(); + const auto& jacobian = linearizer.jacobian(); + linearSolver_.prepare(jacobian, residual); + linearSolver_.getResidual(residual); + solveTimer_.stop(); + + // The preSolve_() method usually computes the errors, but it can do something else in addition. + // TODO: should its costs be counted to the linearization or to the update? + updateTimer_.start(); + preSolve_(currentSolution, residual); + updateTimer_.stop(); + + // Check convergence criteria + if (converged()) { + // Tell the implementation that we're done with this iteration + prePostProcessTimer_.start(); + endIteration_(); + prePostProcessTimer_.stop(); + + break; + } + + // Solve A x = b, where b is the residual, A is its Jacobian and x is the update of the solution + solveTimer_.start(); + solutionUpdate = 0.0; + const bool conv = linearSolver_.solve(solutionUpdate); + solveTimer_.stop(); + + if (!conv) { + solveTimer_.stop(); + if (verbosity_() > 0) { + std::cout << "TPSA: Linear solver did not converge!" << std::endl; + } + + prePostProcessTimer_.start(); + failed_(); + prePostProcessTimer_.stop(); + + return false; + } + + // Update the current solution with the delta + updateTimer_.start(); + update_(nextSolution, currentSolution, solutionUpdate, residual); + updateTimer_.stop(); + + // End of iteration calculations + prePostProcessTimer_.start(); + endIteration_(); + prePostProcessTimer_.stop(); + } + } + catch (const Dune::Exception& e) + { + if (verbosity_() > 0) { + std::cout << "TPSA: Newton method caught exception: \"" + << e.what() << "\"\n" << std::flush; + } + + prePostProcessTimer_.start(); + failed_(); + prePostProcessTimer_.stop(); + + return false; + } + catch (const NumericalProblem& e) + { + if (verbosity_() > 0) { + std::cout << "TPSA: Newton method caught exception: \"" + << e.what() << "\"\n" << std::flush; + } + + prePostProcessTimer_.start(); + failed_(); + prePostProcessTimer_.stop(); + + return false; + } + + // print the timing summary of the time step + if (verbosity_() > 0) { + Scalar elapsedTot = + linearizeTimer_.realTimeElapsed() + + solveTimer_.realTimeElapsed() + + updateTimer_.realTimeElapsed(); + const auto default_precision{std::cout.precision()}; + std::cout << std::setprecision(2) + << "TPSA: " + << "Newton iter = " << numIterations() << " (error=" + << error_ << ") | " + << "linearization = " + << linearizeTimer_.realTimeElapsed() << "s (" + << 100 * linearizeTimer_.realTimeElapsed() / elapsedTot << "%) | " + << "solve = " + << solveTimer_.realTimeElapsed() << "s (" + << 100 * solveTimer_.realTimeElapsed() / elapsedTot << "%) | " + << "update = " + << updateTimer_.realTimeElapsed() << "s (" + << 100 * updateTimer_.realTimeElapsed() / elapsedTot << "%)" + << "\n" << std::flush; + std::cout << std::setprecision(default_precision); // restore default output width + } + + // if we're not converged, tell the implementation that we've failed + if (!converged()) { + prePostProcessTimer_.start(); + failed_(); + if (verbosity_() > 0) { + std::cout << "TPSA: Newton iterations did not converge!" << std::endl; + } + prePostProcessTimer_.stop(); + return false; + } + return true; + } + + /*! + * \brief Returns true if the error of the solution is below the tolerance. + * + * \returns Bool indicating if convergence has been achived + */ + bool converged() const + { return error_ <= tolerance(); } + + /*! + * \brief Returns a reference to the object describing the current physical problem. + * + * \returns Reference to problem object + */ + Problem& problem() + { return simulator_.problem(); } + + /*! + * \brief Returns a reference to the object describing the current physical problem. + * + * \returns Reference to problem object + */ + const Problem& problem() const + { return simulator_.problem(); } + + /*! + * \brief Returns a reference to the geomechanics model. + * + * \returns Reference to geomechanics model object + */ + Model& model() + { return simulator_.problem().geoMechModel(); } + + /*! + * \brief Returns a reference to the geomechanics model. + * + * \returns Reference to geomechanics model object + */ + const Model& model() const + { return simulator_.problem().geoMechModel(); } + + /*! + * \brief Returns the linear solver backend object for external use. + * + * \returns Reference to linear solver object + */ + LinearSolverBackend& linearSolver() + { return linearSolver_; } + + /*! + * \brief Returns the linear solver backend object for external use. + * + * \returns Reference to linear solver object + */ + const LinearSolverBackend& linearSolver() const + { return linearSolver_; } + + /*! + * \brief Returns the number of iterations done since the Newton method was invoked. + * + * \returns Number of Newton iteratinos + */ + int numIterations() const + { return numIterations_; } + + /*! + * \brief Returns the number of linearizations that has done since the Newton method was invoked. + * + * \returns Number of linearizations + */ + int numLinearizations() const + { return numLinearizations_; } + + /*! + * \brief Return the current tolerance at which the Newton method considers itself to be converged. + * + * \returns Tolerance for Newton error + */ + Scalar tolerance() const + { return params_.tolerance_; } + + /*! + * \brief Returns minimum number of Newton iterations used + * + * \returns Minimum number of Newton iterations + */ + Scalar minIterations() const + { return params_.minIterations_; } + + /*! + * \brief Return post-process timer + * + * \returns Reference to post-process timer object + */ + const Timer& prePostProcessTimer() const + { return prePostProcessTimer_; } + + /*! + * \brief Return linearization timer + * + * \returns Reference to linearization timer object + */ + const Timer& linearizeTimer() const + { return linearizeTimer_; } + + /*! + * \brief Return linear solver timer + * + * \returns Reference to linear solver timer object + */ + const Timer& solveTimer() const + { return solveTimer_; } + + /*! + * \brief Return solution update timer + * + * \returns Reference to solution update timer object + */ + const Timer& updateTimer() const + { return updateTimer_; } + +protected: + /*! + * \brief Verbosity level of Newton print messages + * + * \returns Level indicating verbosity + * + * Newton procedure verbosity: 0=none, 1=basic, 2=all + */ + int verbosity_() const + { return simulator_.gridView().comm().rank() == 0 ? params_.verbosity_ : 0; } + + /*! + * \brief Called before the Newton method is applied to an non-linear system of equations. + */ + void begin_() + { + numIterations_ = 0; + numLinearizations_ = 0; + error_ = 1e100; + } + + /*! + * \brief Calculations at the beginning of a Newton iteration + */ + void beginIteration_() + { + // Reset last Newton error to previous + lastError_ = error_; + } + + /*! + * \brief Linearize the global non-linear system of equations associated with the spatial domain. + */ + void linearizeDomain_() + { + model().linearizer().linearizeDomain(); + ++numLinearizations_; + } + + /*! + * \brief Compute error before a Newton iteration + * + * \param currentSolution Current solution vector + * \param currentResidual Current residual vector + */ + void preSolve_(const SolutionVector& /*currentSolution*/, + const GlobalEqVector& currentResidual) + { + lastError_ = error_; + Scalar newtonMaxError = params_.maxError_; + + // Calculate the error as the maximum weighted tolerance of the solution's residual + error_ = 0; + const auto& elemMapper = simulator_.model().elementMapper(); + for (const auto& elem : elements(simulator_.gridView(), Dune::Partitions::interior)) { + unsigned dofIdx = elemMapper.index(elem); + const auto& r = currentResidual[dofIdx]; + for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) { + error_ = max(std::abs(r[eqIdx] * model().eqWeight(dofIdx, eqIdx)), error_); + } + } + + // Take the other processes into account + error_ = simulator_.gridView().comm().max(error_); + + // Make sure that the error never grows beyond the maximum allowed one + if (error_ > newtonMaxError) { + throw NumericalProblem("TPSA: Newton error " + std::to_string(double(error_)) + + " is larger than maximum allowed error of " + + std::to_string(double(newtonMaxError))); + } + } + + /*! + * \brief Update the current solution with a delta vector. + * + * \param nextSolution The solution vector after the current iteration + * \param currentSolution The solution vector after the last iteration + * \param solutionUpdate The delta vector as calculated by solving the linear system of equations + * \param currentResidual The residual vector of the current Newton-Raphson iteraton + * + * Different update strategies, such as chopped updates can be implemented by overriding this method. The default + * behavior is use the standard Newton-Raphson update strategy, i.e. u^{n+1} = u^n - \Delta u^n + */ + void update_(SolutionVector& nextSolution, + const SolutionVector& currentSolution, + const GlobalEqVector& solutionUpdate, + const GlobalEqVector& currentResidual) + { + // make sure not to swallow non-finite values at this point + if (!std::isfinite(solutionUpdate.one_norm())) { + throw NumericalProblem("TPSA: Non-finite update in Newton!"); + } + + std::size_t numGridDof = model().numGridDof(); + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) { + updatePrimaryVariables_(dofIdx, + nextSolution[dofIdx], + currentSolution[dofIdx], + solutionUpdate[dofIdx], + currentResidual[dofIdx]); + } + + // Update material state + model().updateMaterialState(/*timeIdx=*/0); + } + + /*! + * \brief Update a single primary variables object. + * + * \param nextValue The solution vector after the current iteration + * \param currentValue The solution vector after the last iteration + * \param update The delta vector as calculated by solving the linear system of equations + * \param currentResidual The residual vector of the current Newton-Raphson iteraton + */ + void updatePrimaryVariables_(unsigned /*dofIdx*/, + PrimaryVariables& nextValue, + const PrimaryVariables& currentValue, + const EqVector& update, + const EqVector& /*currentResidual*/) + { + nextValue = currentValue; + nextValue -= update; + } + + /*! + * \brief Indicates that one Newton iteration was finished. + */ + void endIteration_() + { + // Increase Newton iterations + ++numIterations_; + + // Output error info + if (verbosity_() > 1) { + std::cout << "TPSA: End Newton iteration " << numIterations_ << "" + << " with error = " << error_ + << std::endl; + } + } + + /*! + * \brief Returns true iff another Newton iteration should be done. + * + * \returns Bool indicating if Newton iterations should continue + */ + bool proceed_() const + { + if (numIterations() < params_.minIterations_) { + return true; + } + else if (converged()) { + // we are below the specified tolerance, so we don't have to + // do more iterations + return false; + } + else if (numIterations() >= params_.maxIterations_) { + // we have exceeded the allowed number of steps. If the + // error was reduced by a factor of at least 4, + // in the last iterations we proceed even if we are above + // the maximum number of steps + return error_ * 4.0 < lastError_; + } + + return true; + } + + /*! + * \brief Called if the Newton method broke down. + */ + void failed_() + { numIterations_ = params_.targetIterations_ * 2; } + + Simulator& simulator_; + LinearSolverBackend linearSolver_; + + Timer prePostProcessTimer_; + Timer linearizeTimer_; + Timer solveTimer_; + Timer updateTimer_; + + Scalar error_; + Scalar lastError_; + TpsaNewtonMethodParams params_; + + int numIterations_; + int numLinearizations_; +}; // class TpsaNewtonMethod + +} // namespace Opm + +#endif diff --git a/opm/models/tpsa/tpsanewtonmethodparams.cpp b/opm/models/tpsa/tpsanewtonmethodparams.cpp new file mode 100644 index 0000000000..789d796d60 --- /dev/null +++ b/opm/models/tpsa/tpsanewtonmethodparams.cpp @@ -0,0 +1,79 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include "config.h" + +#include +#include + +#if HAVE_QUAD +#include +#endif + + +namespace Opm { + +/*! +* \brief Register runtime parameters for TPSA Newton method +*/ +template +void TpsaNewtonMethodParams::registerParameters() +{ + Parameters::Register + ("Verbosity level of TPSA Newton solver: 0 (=none), 1 (=basic), 2 (=all)"); + Parameters::Register + ("The 'optimum' number of TPSA Newton iterations"); + Parameters::Register + ("The maximum number of TPSA Newton iterations"); + Parameters::Register + ("The minimum number of TPSA Newton iterations"); + Parameters::Register> + ("The maximum raw error tolerated by the TPSA Newton method for considering a solution to be converged"); + Parameters::Register> + ("The maximum error tolerated by the TPSA Newton method to which does not cause an abort"); +} + +/*! +* \brief Read and internalize runtime parameters for TPSA Newton method +*/ +template +void TpsaNewtonMethodParams::read() +{ + verbosity_ = Parameters::Get(); + targetIterations_ = Parameters::Get(); + minIterations_ = Parameters::Get(); + maxIterations_ = Parameters::Get(); + tolerance_ = Parameters::Get>(); + maxError_ = Parameters::Get>(); +} + +template struct TpsaNewtonMethodParams; + +#if FLOW_INSTANTIATE_FLOAT +template struct TpsaNewtonMethodParams; +#endif + +#if HAVE_QUAD +template struct TpsaNewtonMethodParams; +#endif + +} // namespace Opm diff --git a/opm/models/tpsa/tpsanewtonmethodparams.hpp b/opm/models/tpsa/tpsanewtonmethodparams.hpp new file mode 100644 index 0000000000..9885b86bd0 --- /dev/null +++ b/opm/models/tpsa/tpsanewtonmethodparams.hpp @@ -0,0 +1,75 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_NEWTON_METHOD_PARAMS_HPP +#define TPSA_NEWTON_METHOD_PARAMS_HPP + + +namespace Opm::Parameters { + +// The maximum error which may occur in a simulation before the Newton method for the time step is aborted +template +struct TpsaNewtonMaxError { static constexpr Scalar value = 1e100; }; + +// Number of maximum iterations for the Newton method +struct TpsaNewtonMaxIterations { static constexpr int value = 20; }; + +// Number of minimum iterations for the Newton method +struct TpsaNewtonMinIterations { static constexpr int value = 1; }; + +// Target number of iterations +struct TpsaNewtonTargetIterations { static constexpr int value = 10; }; + +// Convergence tolerance +template +struct TpsaNewtonTolerance { static constexpr Scalar value = 1e-3; }; + +// Specifies verbosity of print messages +struct TpsaNewtonVerbosity { static constexpr int value = 1; }; + +} // namespace Opm::Parameters + +namespace Opm { + +/*! +* \brief Struct holding the parameters for TpsaNewtonMethod. +*/ +template +struct TpsaNewtonMethodParams +{ + static void registerParameters(); + void read(); + + int verbosity_; + bool writeConvergence_; + int targetIterations_; + int minIterations_; + int maxIterations_; + Scalar tolerance_; + Scalar maxError_; +}; // struct TpsaNewtonMethodParams + +} // namespace Opm + +#endif diff --git a/opm/simulators/flow/BlackoilModelTPSA.hpp b/opm/simulators/flow/BlackoilModelTPSA.hpp new file mode 100644 index 0000000000..620761a43a --- /dev/null +++ b/opm/simulators/flow/BlackoilModelTPSA.hpp @@ -0,0 +1,230 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef BLACK_OIL_MODEL_TPSA_HPP +#define BLACK_OIL_MODEL_TPSA_HPP + +#include +#include + +#include + +#include +#include + +#include + + +namespace Opm { + +/*! +* \brief Black oil model for coupling Flow simulations with TPSA geomechanics +*/ +template +class BlackoilModelTPSA : public BlackoilModel +{ + using ParentType = BlackoilModel; + + using Scalar = GetPropType; + using Simulator = GetPropType; + +public: + using ModelParameters = typename ParentType::ModelParameters; + + /*! + * \brief Constructor + * + * \param simulator Reference to simulator object + * \param param Reference to parameters for model + * \param well_model Refenerence to well model + * \param terminal_output Bool for terminal output + */ + explicit BlackoilModelTPSA(Simulator& simulator, + const ModelParameters& param, + BlackoilWellModel& well_model, + const bool terminal_output) + : ParentType(simulator, param, well_model, terminal_output) + {} + + /*! + * \brief Perform a nonlinear iteration updating Flow and TPSA geomechanics + * + * \param timer Simulation timer + * \param nonlinear_solver Nonlinear solver type + * \returns Report for simulator performance + * + * \note Strategies of coupling Flow and TPSA currently implemented: + * \li fixed-stress: fixed-stress algorithm, i.e. iteratively solving Flow and TPSA equations in sequence + * \li lagged: one-way coupling where Flow is solved with TPSA info from previous time step + */ + template + SimulatorReportSingle nonlinearIteration(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + SimulatorReportSingle report {}; + const auto& problem = this->simulator_.problem(); + if (problem.fixedStressScheme()) { + report = nonlinearIterationFixedStressTPSA(timer, nonlinear_solver); + } + else if (problem.laggedScheme()) { + report = nonlinearIterationLaggedTPSA(timer, nonlinear_solver); + } + else { + std::string msg = "Unknown Flow-TPSA coupling scheme!"; + OpmLog::error(msg); + throw std::runtime_error(msg); + } + return report; + } + + /*! + * \brief Perform a nonlinear iteration updating Flow and TPSA geomechanics in a fixed-stress, iterative loop + * + * \param timer Simulation timer + * \param nonlinear_solver Nonlinear solver type + * \returns Report for simulator performance + * + */ + template + SimulatorReportSingle nonlinearIterationFixedStressTPSA(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + // Runtime parameters + const auto& [minSeqIter, maxSeqIter] = this->simulator_.problem().fixedStressParameters(); + SimulatorReportSingle reportFlow; + + // Max. no. of fixed-stress iterations reached: warn and move on + if (seqIter_ >= maxSeqIter) { + // Warning + std::string msg = fmt::format("TPSA: Fixed-stress scheme reached max iterations (={})!", maxSeqIter); + OpmLog::warning(msg); + + // Return true Flow convergence to move to next time step and reset other variables + reportFlow.converged = true; + seqIter_ = 0; + + return reportFlow; + } + + // Prepare before first iteration + if (seqIter_ == 0) { + this->simulator_.problem().geoMechModel().prepareTPSA(); + } + + // Run Flow nonlinear iteration + reportFlow = ParentType::nonlinearIteration(timer, nonlinear_solver); + + // Solve TPSA equations if: + // (i) Flow has converged and run more than min number of Newton iterations + // (ii) we have run at least min. number of fixed-stress iterations + const auto iteration = this->simulator_.problem().iterationContext().iteration(); + if (reportFlow.converged && iteration >= this->param_.newton_min_iter_) { + // Solve TPSA equations + bool tpsaConv = solveTpsaEquations(); + ++seqIter_; + + // Fixed-stress convergence check: + // If the initial residual error, hence check for no. linearizations == 1, was small enough, we have + // convergence in the fixed-stress iterations + if (tpsaConv + && this->simulator_.problem().geoMechModel().newtonMethod().numLinearizations() == 1 + && seqIter_ >= minSeqIter) { + // Info + std::string msg = fmt::format("TPSA: Fixed-stress scheme converged in {} iterations", seqIter_); + OpmLog::info(msg); + + // Reset + seqIter_ = 0; + + return reportFlow; + } + // Throw error if TPSA did not converge. Will force time step cuts in the outer Flow loop. + else if (!tpsaConv) { + // Reset + seqIter_ = 0; + + throw std::runtime_error("TPSA: Fixed stress scheme update failed!"); + } + + // Return Flow convergence false to do another fixed-stress iteration + reportFlow.converged = false; + } + + return reportFlow; + } + + /*! + * \brief Perform a nonlinear iteration updating Flow and TPSA geomechanics in a lagged scheme + * + * \param iteration Flow nonlinear iteration + * \param timer Simulation timer + * \param nonlinear_solver Nonlinear solver type + * \returns Report for simulator performance + * + */ + template + SimulatorReportSingle nonlinearIterationLaggedTPSA(const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + // Run Flow nonlinear iteration + auto reportFlow = ParentType::nonlinearIteration(timer, nonlinear_solver); + + // Update TPSA geomechanics from successful Flow run + const auto iteration = this->simulator_.problem().iterationContext().iteration(); + if (reportFlow.converged && iteration >= this->param_.newton_min_iter_) { + // Prepare before TPSA solve + this->simulator_.problem().geoMechModel().prepareTPSA(); + + // Solve TPSA equations + bool tpsaConv = solveTpsaEquations(); + + // Throw error if TPSA did not converge. Will force time step cuts in the outer Flow loop. + if (!tpsaConv) { + throw std::runtime_error("TPSA: Lagged scheme update failed!"); + } + } + + return reportFlow; + } + + /*! + * \brief Solve TPSA geomechanics equations + * + * \returns Bool indicating TPSA convergence + * + * \note Calls Newton method for TPSA + */ + bool solveTpsaEquations() + { + // Run Newthon method for TPSA equations + return this->simulator_.problem().geoMechModel().newtonMethod().apply(); + } + +private: + int seqIter_{0}; +}; // class BlackoilModelTPSA + +} // namespace Opm + +#endif diff --git a/opm/simulators/flow/FacePropertiesTPSA.cpp b/opm/simulators/flow/FacePropertiesTPSA.cpp new file mode 100644 index 0000000000..d2356cf5a7 --- /dev/null +++ b/opm/simulators/flow/FacePropertiesTPSA.cpp @@ -0,0 +1,76 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include + +#include + +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif + + +namespace Opm { + +#define INSTANTIATE_TYPE(T) \ + template class FacePropertiesTPSA>, \ + Dune::MultipleCodimMultipleGeomTypeMapper< \ + Dune::GridView< \ + Dune::DefaultLeafGridViewTraits>>, \ + Dune::CartesianIndexMapper, \ + T>; + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +#ifdef HAVE_DUNE_FEM +using GV = Dune::Fem::AdaptiveLeafGridPart; +template class FacePropertiesTPSA, + Dune::CartesianIndexMapper, + double>; + +#if FLOW_INSTANTIATE_FLOAT +template class FacePropertiesTPSA, + Dune::CartesianIndexMapper, + float>; +#endif + +#endif // HAVE_DUNE_FEM + +} // namespace Opm diff --git a/opm/simulators/flow/FacePropertiesTPSA.hpp b/opm/simulators/flow/FacePropertiesTPSA.hpp new file mode 100644 index 0000000000..fb71867593 --- /dev/null +++ b/opm/simulators/flow/FacePropertiesTPSA.hpp @@ -0,0 +1,141 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_FACE_PROPERTIES_HPP +#define TPSA_FACE_PROPERTIES_HPP + +#include + +#include + +#include +#include +#include +#include + + +namespace Opm { + +class EclipseState; + +/*! +* \brief Cell face properties needed in TPSA equation calculations +* +* Similar calculations as done in Transmissibility class for TPFA +*/ +template +class FacePropertiesTPSA { + + enum { dimWorld = GridView::dimensionworld }; + +public: + using DimVector = Dune::FieldVector; + + FacePropertiesTPSA(const EclipseState& eclState, + const GridView& gridView, + const CartesianIndexMapper& cartMapper, + const Grid& grid, + std::function(int)> centroids); + + void finishInit(); + + void update(); + + Scalar weightAverage(unsigned elemIdx1, unsigned elemIdx2) const; + Scalar weightAverageBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const; + Scalar weightProduct(unsigned elemIdx1, unsigned elemIdx2) const; + Scalar weightProductBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const; + Scalar normalDistance(unsigned elemIdx1, unsigned elemIdx2) const; + Scalar normalDistanceBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const; + DimVector cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2); + const DimVector& cellFaceNormalBoundary(unsigned elemIdx1, unsigned boundaryFaceIdx) const; + + /*! + * \brief Return shear modulus of an element + */ + const Scalar shearModulus(unsigned elemIdx) const + { return sModulus_[elemIdx]; } + + +protected: + struct FaceInfo + { + DimVector faceCenter; + int faceIdx; + unsigned elemIdx; + unsigned cartElemIdx; + }; + + template + void computeCellProperties(const Intersection& intersection, + FaceInfo& inside, + FaceInfo& outside, + DimVector& faceNormal, + /*isCpGrid=*/std::false_type) const; + + template + void computeCellProperties(const Intersection& intersection, + FaceInfo& inside, + FaceInfo& outside, + DimVector& faceNormal, + /*isCpGrid=*/std::true_type) const; + + Scalar computeDistance_(const DimVector& distVec, const DimVector& faceNormal); + + Scalar computeWeight_(const Scalar distance, const Scalar smod); + + DimVector distanceVector_(const DimVector& faceCenter, const unsigned& cellIdx) const; + + void extractSModulus_(); + + std::vector sModulus_; + std::unordered_map weightsAvg_; + std::unordered_map weightsProd_; + std::unordered_map distance_; + std::unordered_map faceNormal_; + + std::map, Scalar> weightsAvgBoundary_; + std::map, Scalar> weightsProdBoundary_; + std::map, Scalar> distanceBoundary_; + std::map, DimVector> faceNormalBoundary_; + + const EclipseState& eclState_; + const GridView& gridView_; + const CartesianIndexMapper& cartMapper_; + const Grid& grid_; + std::function(int)> centroids_; + std::vector> centroids_cache_; + + const LookUpData lookUpData_; + const LookUpCartesianData lookUpCartesianData_; + +}; // FacePropertiesTPSA + +namespace details { + std::uint64_t isIdTPSA(std::uint32_t elemIdx1, std::uint32_t elemIdx2); +} // namespace details + +} // namespace Opm + +#endif diff --git a/opm/simulators/flow/FacePropertiesTPSA_impl.hpp b/opm/simulators/flow/FacePropertiesTPSA_impl.hpp new file mode 100644 index 0000000000..1d3dacf4b7 --- /dev/null +++ b/opm/simulators/flow/FacePropertiesTPSA_impl.hpp @@ -0,0 +1,604 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_FACE_PROPERTIES_IMPL_HPP +#define TPSA_FACE_PROPERTIES_IMPL_HPP + +#ifndef TPSA_FACE_PROPERTIES_HPP +#include +#include +#endif + +#include + +#include + +#include + +#include + +#include + +#include +#include +#include + + +namespace Opm { + +// Copied from Opm::Transmissibility class +namespace details { + constexpr unsigned elemIdxShift = 32; // bits + + std::uint64_t isIdTPSA(std::uint32_t elemIdx1, std::uint32_t elemIdx2) + { + const std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2); + const std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2); + + return (elemBIdx << elemIdxShift) + elemAIdx; + } +} // namespace Opm::details + +// ///// +// Public functions +// //// +/*! +* \brief Constructor +* +* \param eclState Eclipse state made from a deck +* \param gridView The view on the DUNE grid which ought to be used (normally the leaf grid view) +* \param cartMapper The cartesian index mapper for lookup of cartesian indices +* \param grid The grid to lookup cell properties +* \param centroids Function to lookup cell centroids +*/ +template +FacePropertiesTPSA:: +FacePropertiesTPSA(const EclipseState& eclState, + const GridView& gridView, + const CartesianIndexMapper& cartMapper, + const Grid& grid, + std::function(int)> centroids) + : eclState_(eclState) + , gridView_(gridView) + , cartMapper_(cartMapper) + , grid_(grid) + , centroids_(centroids) + , lookUpData_(gridView) + , lookUpCartesianData_(gridView, cartMapper) +{ } + +/*! +* \brief Compute TPSA face properties +*/ +template +void FacePropertiesTPSA:: +finishInit() +{ + update(); +} + +/*! +* \brief Compute TPSA face properties +*/ +template +void FacePropertiesTPSA:: +update() +{ + // Number of elements + ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); + unsigned numElements = elemMapper.size(); + + // Extract shear modulus (Lame's sec. params) + extractSModulus_(); + + // Init. containers + // Note (from Transmissibility::update): Reserving some space in the hashmap upfront saves quite a bit of time + // because resizes are costly for hashmaps and there would be quite a few of them if we would not have a rough idea + // of how large the final map will be (the rough idea is a conforming Cartesian grid). + const int num_threads = ThreadManager::maxThreads(); + weightsAvg_.clear(); + weightsProd_.clear(); + distance_.clear(); + faceNormal_.clear(); + if (num_threads == 1) { + weightsAvg_.reserve(numElements * 3 * 1.05); + weightsProd_.reserve(numElements * 3 * 1.05); + distance_.reserve(numElements * 3 * 1.05); + faceNormal_.reserve(numElements * 3 * 1.05); + } + weightsAvgBoundary_.clear(); + weightsProdBoundary_.clear(); + distanceBoundary_.clear(); + faceNormalBoundary_.clear(); + + // Fill the centroids cache to avoid repeated calculations in loops below + centroids_cache_.resize(gridView_.size(0)); + for (const auto& elem : elements(gridView_)) { + const unsigned elemIdx = elemMapper.index(elem); + centroids_cache_[elemIdx] = centroids_(elemIdx); + } + + // Initialize thread safe insert_or_assign for face properties in the grid and separate for boundaries + ThreadSafeMapBuilder weightsAvgMap(weightsAvg_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder weightsProdMap(weightsProd_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder distanceMap(distance_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder faceNormalMap(faceNormal_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + + ThreadSafeMapBuilder weightsAvgBoundaryMap(weightsAvgBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder weightsProdBoundaryMap(weightsProdBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder distanceBoundaryMap(distanceBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder faceNormalBoundaryMap(faceNormalBoundary_, num_threads, MapBuilderInsertionMode::Insert_Or_Assign); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + // Loop over grid element an compute face properties + for (const auto& chunk : ElementChunks(gridView_, Dune::Partitions::all, num_threads)) { + for (const auto& elem : chunk) { + // Init. face info for inside/outside cells + FaceInfo inside; + FaceInfo outside; + DimVector faceNormal; + + // Set inside info + inside.elemIdx = elemMapper.index(elem); + inside.cartElemIdx = lookUpCartesianData_. + template getFieldPropCartesianIdx(inside.elemIdx); + + // Loop over intersection to neighboring cells + unsigned boundaryIsIdx = 0; + for (const auto& intersection : intersections(gridView_, elem)) { + // Handle grid boundary + if (intersection.boundary()) { + // One-sided cell properties + const auto& geometry = intersection.geometry(); + inside.faceCenter = geometry.center(); + faceNormal = intersection.centerUnitOuterNormal(); + + // Face properties on boundary + const auto index_pair = std::make_pair(inside.elemIdx, boundaryIsIdx); + Scalar distBound = computeDistance_(distanceVector_(inside.faceCenter, inside.elemIdx), faceNormal); + distanceBoundaryMap.insert_or_assign(index_pair, distBound); + + Scalar weightsAvgBound = 1.0; // w_j = 0 -> w_avg = wi / (wi + 0) + Scalar weightsProdBound = 0.0; // w_j = 0 -> w_prod = wi * 0 + weightsAvgBoundaryMap.insert_or_assign(index_pair, weightsAvgBound); + weightsProdBoundaryMap.insert_or_assign(index_pair, weightsProdBound); + + faceNormalBoundaryMap.insert_or_assign(index_pair, faceNormal); + + ++boundaryIsIdx; + continue; + } + + // Handle intersection on process boundary (i.e., neighbor on different rank) + if (!intersection.neighbor()) { + ++boundaryIsIdx; + continue; + } + + // Set outside info + const auto& outsideElem = intersection.outside(); + outside.elemIdx = elemMapper.index(outsideElem); + outside.cartElemIdx = lookUpCartesianData_. + template getFieldPropCartesianIdx(outside.elemIdx); + + // Skip intersections that have already been processed + if (std::tie(inside.cartElemIdx, inside.elemIdx) > std::tie(outside.cartElemIdx, outside.elemIdx)) { + continue; + } + + // Face indices for this intersection + inside.faceIdx = intersection.indexInInside(); + outside.faceIdx = intersection.indexInOutside(); + + // Set NNC face properties to zero + if (inside.faceIdx == -1) { + const auto id = details::isIdTPSA(inside.elemIdx, outside.elemIdx); + weightsAvgMap.insert_or_assign(id, 0.0); + distanceMap.insert_or_assign(id, 0.0); + faceNormalMap.insert_or_assign(id, DimVector{0.0, 0.0, 0.0}); + + continue; + } + + // Compute cell properties from grid + typename std::is_same::type isCpGrid; + computeCellProperties(intersection, + inside, + outside, + faceNormal, + isCpGrid); + + // Compute face properties + const auto id = details::isIdTPSA(inside.elemIdx, outside.elemIdx); + Scalar dist_in = computeDistance_(distanceVector_(inside.faceCenter, inside.elemIdx), faceNormal); + Scalar dist_out = computeDistance_(distanceVector_(outside.faceCenter, outside.elemIdx), faceNormal); + distanceMap.insert_or_assign(id, dist_in + dist_out); + + Scalar weight_in = computeWeight_(dist_in, sModulus_[inside.elemIdx]); + Scalar weight_out = computeWeight_(dist_out, sModulus_[outside.elemIdx]); + Scalar weightsAvg = weight_in / (weight_in + weight_out); + Scalar weightsProd = weight_in * weight_out; + weightsAvgMap.insert_or_assign(id, weightsAvg); + weightsProdMap.insert_or_assign(id, weightsProd); + + faceNormalMap.insert_or_assign(id, faceNormal); + } + } + } + + // Clear centroid cache + centroids_cache_.clear(); + +#ifdef _OPENMP +#pragma omp parallel sections +#endif + { +#ifdef _OPENMP +#pragma omp section +#endif + weightsAvgMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + weightsProdMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + distanceMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + faceNormalMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + weightsAvgBoundaryMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + weightsProdBoundaryMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + distanceBoundaryMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + faceNormalBoundaryMap.finalize(); + } +} + +/*! +* \brief Average (half-)weight at interface between two elements +* +* \param elemIdx1 Cell index 1 +* \param elemIdx2 Cell index 2 +* \returns Average weight +*/ +template +Scalar FacePropertiesTPSA:: +weightAverage(unsigned elemIdx1, unsigned elemIdx2) const +{ + auto tmp_whgt = weightsAvg_.at(details::isIdTPSA(elemIdx1, elemIdx2)); + + auto cartIdx1 = lookUpCartesianData_. + template getFieldPropCartesianIdx(elemIdx1); + auto cartIdx2 = lookUpCartesianData_. + template getFieldPropCartesianIdx(elemIdx2); + if (cartIdx1 < cartIdx2) { + return tmp_whgt; + } + else { + return 1.0 - tmp_whgt; + } +} + +/*! +* \brief Average (half-)weight at boundary interface +* +* \param elemIdx Cell index +* \param boundaryFaceIdx Face index +* \returns Average weight at boundary +*/ +template +Scalar FacePropertiesTPSA:: +weightAverageBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const +{ + return weightsAvgBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); +} + +/*! +* \brief Product of weights at interface between two elements +* +* \param elemIdx1 Cell index 1 +* \param elemIdx2 Cell index 2 +* \returns Weight product +*/ +template +Scalar FacePropertiesTPSA:: +weightProduct(unsigned elemIdx1, unsigned elemIdx2) const +{ + return weightsProd_.at(details::isIdTPSA(elemIdx1, elemIdx2)); +} + +/*! +* \brief Product of weights at boundary interface +* +* \param elemIdx Cell index +* \param boundaryFaceIdx Face index +* \returns Weight product at boundary +*/ +template +Scalar FacePropertiesTPSA:: +weightProductBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const +{ + return weightsProdBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); +} + +/*! +* \brief Distance between two elements +* +* \param elemIdx1 Cell index 1 +* \param elemIdx2 Cell index 2 +* \returns Distance +*/ +template +Scalar FacePropertiesTPSA:: +normalDistance(unsigned elemIdx1, unsigned elemIdx2) const +{ + return distance_.at(details::isIdTPSA(elemIdx1, elemIdx2)); +} + +/*! +* \brief Distance to boundary interface +* +* \param elemIdx Cell index +* \param boundaryFaceIdx Face index +* \returns Distance to boundary +*/ +template +Scalar FacePropertiesTPSA:: +normalDistanceBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const +{ + return distanceBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); +} + +/*! +* \brief Cell face normal at interface between two elements +* +* \param elemIdx1 Cell index 1 +* \param elemIdx2 Cell index 2 +* \returns Face normals +* +* \note It is assumed that normals point from lower index to higher index! +*/ +template +typename FacePropertiesTPSA::DimVector +FacePropertiesTPSA:: +cellFaceNormal(unsigned elemIdx1, unsigned elemIdx2) +{ + auto cartIdx1 = lookUpCartesianData_. + template getFieldPropCartesianIdx(elemIdx1); + auto cartIdx2 = lookUpCartesianData_. + template getFieldPropCartesianIdx(elemIdx2); + int sign = (cartIdx1 < cartIdx2) ? 1 : -1; + return sign * faceNormal_.at(details::isIdTPSA(elemIdx1, elemIdx2)); +} + +/*! +* \brief Cell face normal of boundary interface +* +* \param elemIdx Cell index +* \param boundaryFaceIdx Face index +* \returns Face normals +* +* \note Boundary normals points outwards +*/ +template +const typename FacePropertiesTPSA::DimVector& +FacePropertiesTPSA:: +cellFaceNormalBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const +{ + return faceNormalBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); +} + +// ///// +// Protected functions +// //// +/*! +* \brief Compute normal distance from cell center to face center +* +* \param distVec Distance vector from cell center to face center +* \param faceNormal Face (unit) normal vector +* \returns Distance cell center -> face center +*/ +template +Scalar FacePropertiesTPSA:: +computeDistance_(const DimVector& distVec, const DimVector& faceNormal) +{ + return std::abs(Dune::dot(faceNormal, distVec)); +} + +/*! +* \brief Compute face properties from general DUNE grid +* +* \param intersection Info on cell intersection +* \param inside Info describing inside face +* \param outside Info describing outside face +* \param faceNormal Face (unit) normal vector +* +* \warning Not implemented! +*/ +template +template +void FacePropertiesTPSA:: +computeCellProperties(const Intersection& intersection, + FaceInfo& inside, + FaceInfo& outside, + DimVector& faceNormal, + /*isCpGrid=*/std::false_type) const +{ + throw std::runtime_error("TPSA not implemented for DUNE grid types other than CpGrid!"); +} + +/*! +* \brief Compute face properties from DUNE CpGrid +* +* \param intersection Info on cell intersection +* \param inside Info describing inside face +* \param outside Info describing outside face +* \param faceNormal Face (unit) normal +* +* \warning LGR computations not implemented! +*/ +template +template +void FacePropertiesTPSA:: +computeCellProperties(const Intersection& intersection, + FaceInfo& inside, + FaceInfo& outside, + DimVector& faceNormal, + /*isCpGrid=*/std::true_type) const +{ + int faceIdx = intersection.id(); + if (grid_.maxLevel() == 0) { + // Face center coordinates + inside.faceCenter = grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection); + outside.faceCenter = grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection); + + // Face normal, ensuring it points from cell with lower to higher (global grid) index + faceNormal = grid_.faceNormal(faceIdx); + auto cartFaceCell0 = lookUpCartesianData_. + template getFieldPropCartesianIdx(grid_.faceCell(faceIdx, 0)); + auto cartFaceCell1 = lookUpCartesianData_. + template getFieldPropCartesianIdx(grid_.faceCell(faceIdx, 1)); + if (cartFaceCell0 > cartFaceCell1) { + faceNormal *= -1; + } + } + else { + throw std::runtime_error("TPSA not implemented with LGR"); + } + +} + +/*! +* \brief Compute weight ratio between distance and shear modulus +* +* \param distance Normal distance from cell to face center +* \param smod Shear modulus +* \returns Weight ratio +*/ +template +Scalar FacePropertiesTPSA:: +computeWeight_(const Scalar distance, const Scalar smod) +{ + return distance / smod; +} + +/*! +* \brief Distance vector from cell center to face center +* +* \param faceCenter Face center coordinates +* \param cellIdx Cell index +* \returns Distance vector cell center -> face center +*/ +template +typename FacePropertiesTPSA::DimVector +FacePropertiesTPSA:: +distanceVector_(const DimVector& faceCenter, const unsigned& cellIdx) const +{ + const auto& cellCenter = centroids_cache_.empty() ? centroids_(cellIdx) + : centroids_cache_[cellIdx]; + DimVector x = faceCenter; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) { + x[dimIdx] -= cellCenter[dimIdx]; + } + + return x; +} + +/*! +* \brief Extract shear modulus from eclState +* +* \note (from Transmissibility::extractPorosity()): +* All arrays provided by eclState are one-per-cell of "uncompressed" grid, whereas the simulation grid might remove a +* few elements (e.g. because it is distributed over several processes). +*/ +template +void FacePropertiesTPSA:: +extractSModulus_() +{ + unsigned numElem = gridView_.size(/*codim=*/0); + sModulus_.resize(numElem); + + const auto& fp = eclState_.fieldProps(); + std::vector sModulusData; + if (fp.has_double("SMODULUS")) { + sModulusData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "SMODULUS"); + } + else if (fp.has_double("YMODULE") && fp.has_double("LAME")) { + // Convert from Young's modulus and Lame's first parameter + const std::vector& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE"); + const std::vector& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME"); + sModulusData.resize(numElem); + for (std::size_t i = 0; i < sModulusData.size(); ++i) { + const double r = std::sqrt(ymodulus[i] * ymodulus[i] + 9 * lameParam[i] * lameParam[i] + + 2 * ymodulus[i] * lameParam[i]); + sModulusData[i] = (ymodulus[i] - 3 * lameParam[i] + r) / 4.0; + } + } + else if (fp.has_double("YMODULE") && fp.has_double("PRATIO")) { + const std::vector& ymodulus = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "YMODULE"); + const std::vector& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO"); + sModulusData.resize(numElem); + for (std::size_t i = 0; i < sModulusData.size(); ++i) { + sModulusData[i] = ymodulus[i] / (2 * (1 + pratio[i])); + } + } + else if (fp.has_double("LAME") && fp.has_double("PRATIO")) { + const std::vector& lameParam = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "LAME"); + const std::vector& pratio = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PRATIO"); + sModulusData.resize(numElem); + for (std::size_t i = 0; i < sModulusData.size(); ++i) { + sModulusData[i] = lameParam[i] * (1 - 2 * pratio[i]) / (2 * pratio[i]); + } + } + else { + throw std::logic_error("Cannot read shear modulus data from ecl state, SMODULUS keyword missing, " + "and one of the following keyword pairs are missing for conversion: " + "(YMODULE, LAME), (YMODULE, PRATIO) and (LAME, PRATIO)!"); + } + + // Assign shear modulus + for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { + sModulus_[dofIdx] = sModulusData[dofIdx]; + } +} + +} // namespace Opm + +#endif diff --git a/opm/simulators/flow/FlowGenericProblem.hpp b/opm/simulators/flow/FlowGenericProblem.hpp index a739545066..e04c2dc6e7 100644 --- a/opm/simulators/flow/FlowGenericProblem.hpp +++ b/opm/simulators/flow/FlowGenericProblem.hpp @@ -150,6 +150,21 @@ public: return rockFraction_[timeIdx][elementIdx]; }; + /*! + * \brief Returns the rock compressibility of an element due to poroelasticity + */ + Scalar rockBiotComp(unsigned elementIdx) const; + + /*! + * \brief Direct access to Lame's second parameter in an element + */ + Scalar lame(unsigned elementIdx) const; + + /*! + * \brief Direct access to Biot coefficient in an element + */ + Scalar biotCoeff(unsigned elementIdx) const; + /*! * \brief Sets the porosity of an element * diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index b0f35e15ab..e15d5e351c 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -340,6 +340,70 @@ porosity(unsigned globalSpaceIdx, unsigned timeIdx) const return this->referencePorosity_[timeIdx][globalSpaceIdx]; } +template +typename FlowGenericProblem::Scalar +FlowGenericProblem:: +rockBiotComp(unsigned elementIdx) const +{ + // Additional compressibility of the rock due to Biot poroelasticity + auto biot = biotCoeff(elementIdx); + auto lameParam = lame(elementIdx); + return biot * biot / lameParam; +} + +template +typename FlowGenericProblem::Scalar +FlowGenericProblem:: +lame(unsigned elementIdx) const +{ + Scalar lameParam; + const auto& fp = eclState_.fieldProps(); + if (fp.has_double("LAME")) { + lameParam = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "LAME", elementIdx); + } + else if (fp.has_double("YMODULE") && fp.has_double("SMODULUS")) { + const auto& yModulus = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "YMODULE", elementIdx); + const auto& sModulus = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "SMODULUS", elementIdx); + lameParam = sModulus * (yModulus - 2 * sModulus) / (3 * sModulus - yModulus); + } + else if (fp.has_double("YMODULE") && fp.has_double("PRATIO")) { + const auto& yModulus = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "YMODULE", elementIdx); + const auto& pRatio = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "PRATIO", elementIdx); + lameParam = yModulus * pRatio / ((1 + pRatio) * (1 - 2 * pRatio)); + } + else if (fp.has_double("SMODULUS") && fp.has_double("PRATIO")) { + const auto& sModulus = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "SMODULUS", elementIdx); + const auto& pRatio = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "PRATIO", elementIdx); + lameParam = 2 * sModulus * pRatio / (1 - 2 * pRatio); + } + else { + return 0.0; + } + return lameParam; +} + +template +typename FlowGenericProblem::Scalar +FlowGenericProblem:: +biotCoeff(unsigned elementIdx) const +{ + Scalar biotC; + const auto& fp = eclState_.fieldProps(); + if (fp.has_double("BIOTCOEF")) { + biotC = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "BIOTCOEF", elementIdx); + } + else if (fp.has_double("POELCOEF") && fp.has_double("PRATIO")) { + const auto& poelC = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "POELCOEF", elementIdx); + const auto& pRatio = this->lookUpData_.fieldPropDouble(this->eclState_.fieldProps(), "PRATIO", elementIdx); + biotC = poelC * (1 - pRatio) / (1 - 2 * pRatio); + } + else { + return 0.0; + } + return biotC; +} + + template template void FlowGenericProblem:: diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index 38d8f3c303..73ff334ffd 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -173,6 +173,9 @@ public: using BaseType::helpPreamble; using BaseType::shouldWriteOutput; using BaseType::shouldWriteRestartFile; + using BaseType::rockBiotComp; + using BaseType::lame; + using BaseType::biotCoeff; using BaseType::rockCompressibility; using BaseType::porosity; @@ -715,6 +718,36 @@ public: return this->rockCompressibility(globalSpaceIdx); } + /*! + * \copydoc BlackoilProblem::rockBiotComp + */ + template + Scalar rockBiotComp(const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + return this->rockBiotComp(globalSpaceIdx); + } + + /*! + * \copydoc BlackoilProblem::lame + */ + template + Scalar lame(const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + return this->lame(globalSpaceIdx); + } + + /*! + * \copydoc BlackoilProblem::biotCoeff + */ + template + Scalar biotCoeff(const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + return this->biotCoeff(globalSpaceIdx); + } + /*! * \copydoc BlackoilProblem::rockReferencePressure */ diff --git a/opm/simulators/flow/FlowProblemTPSA.hpp b/opm/simulators/flow/FlowProblemTPSA.hpp new file mode 100644 index 0000000000..936279753b --- /dev/null +++ b/opm/simulators/flow/FlowProblemTPSA.hpp @@ -0,0 +1,536 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef FLOW_PROBLEM_TPSA_HPP +#define FLOW_PROBLEM_TPSA_HPP + +#include + +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include + + +namespace Opm { + +/*! +* \brief Problem for Flow-TPSA coupled simulations +*/ +template +class FlowProblemTPSA : public FlowProblemBlackoil +{ +public: + using ParentType = FlowProblemBlackoil; + + using ElementContext = GetPropType; + using ElementMapper = GetPropType; + using Evaluation = GetPropType; + using FluidSystem = GetPropType; + using GeomechModel = GetPropType; + using Grid = GetPropType; + using GridView = GetPropType; + using Indices = GetPropType; + using RateVector = GetPropType; + using Scalar = GetPropType; + using Simulator = GetPropType; + + enum { dimWorld = GridView::dimensionworld }; + enum { enableMech = getPropValue() }; + enum { historySize = getPropValue() }; + enum { numEq = getPropValue() }; + enum { numPhases = FluidSystem::numPhases }; + + enum { contiRotEqIdx = Indices::contiRotEqIdx }; + enum { contiSolidPresEqIdx = Indices::contiSolidPresEqIdx }; + enum { solidPres0Idx = Indices::solidPres0Idx }; + + using CartesianIndexMapper = Dune::CartesianIndexMapper; + using DimVector = Dune::FieldVector; + using FaceProperties = FacePropertiesTPSA; + using InitialMaterialState = MaterialStateTPSA; + using Toolbox = MathToolbox; + + // /// + // Public functions + // /// + /*! + * \brief Constructor + * + * \param simulator Reference to simulator object + */ + FlowProblemTPSA(Simulator& simulator) + : ParentType(simulator) + , faceProps_(simulator.vanguard().eclState(), + simulator.vanguard().gridView(), + simulator.vanguard().cartesianIndexMapper(), + simulator.vanguard().grid(), + simulator.vanguard().cellCentroids()) + , geoMechModel_(simulator) + { + if constexpr(enableMech) { + // Add VTK TPSA to output module + this->model().addOutputModule(std::make_unique>(simulator)); + + // Sanity check + const auto& mechSolver = simulator.vanguard().eclState().runspec().mechSolver(); + if (!mechSolver.tpsa()) { + std::string msg = "Simulator with Tpsa-geomechanics enabled compile time, but deck does not contain " + "TPSA keyword!"; + OpmLog::error(msg); + throw std::runtime_error(msg); + } + } + else { + // Sanity check + const auto& mechSolver = simulator.vanguard().eclState().runspec().mechSolver(); + if (mechSolver.tpsa()) { + std::string msg = "TPSA keyword in deck, but Tpsa-geomechanics disabled compile-time!"; + OpmLog::error(msg); + throw std::runtime_error(msg); + } + } + } + + /*! + * \brief Register runtime parameters + */ + static void registerParameters() + { + // Register parameters for parent class + ParentType::registerParameters(); + + // Geomech model parameters + GeomechModel::registerParameters(); + + // VTK output parameters + VtkTpsaModule::registerParameters(); + } + + /*! + * \brief Initialize the problem + */ + void finishInit() + { + // FlowProblemBlackoil::finishInit() + ParentType::finishInit(); + + // Read initial conditions and set material state + readInitalConditionsTPSA_(); + + // Calculate face properties + faceProps_.finishInit(); + + // Set equation weights + computeAndSetEqWeights_(); + + // Initialize the TPSA model + geoMechModel_.finishInit(); + } + + /*! + * \brief Set initial solution for the problem + * + * This function is a combination FvBaseDiscretization::applyInitialSolution and FlowProblemBlackoil::initial() + */ + void initialSolutionApplied() + { + // Set up initial solution for the Flow model + ParentType::initialSolutionApplied(); + + // Initialize soultions as zero + auto& uCur = geoMechModel_.solution(/*timeIdx=*/0); + uCur = Scalar(0.0); + + // Loop through grid and set initial solution from material state + ElementContext elemCtx(this->simulator()); + for (const auto& elem : elements(this->gridView())) { + // Ignore everything which is not in the interior if the current process' piece of the grid + if (elem.partitionType() != Dune::InteriorEntity) { + continue; + } + + // Loop over sub control volumes and set initial solutions + elemCtx.updateStencil(elem); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + auto& priVars = uCur[globalIdx]; + priVars.assignNaive(initialMaterialState_[globalIdx]); + priVars.checkDefined(); + } + } + + // synchronize the ghost/overlapping DOFs (if necessary) + geoMechModel_.syncOverlap(); + + // Set history solutions to the initial solution. + for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx) { + geoMechModel_.solution(timeIdx) = geoMechModel_.solution(/*timeIdx=*/0); + } + + // Set material state + geoMechModel_.updateMaterialState(/*timeIdx=*/0); + } + + /*! + * \brief Compute weights to rescale the TPSA equations + */ + void computeAndSetEqWeights_() + { + // Average shear modulus over domain + Scalar avgSmodulus = 0.0; + const auto& gridView = this->gridView(); + ElementContext elemCtx(this->simulator()); + for(const auto& elem: elements(gridView, Dune::Partitions::interior)) { + elemCtx.updatePrimaryStencil(elem); + int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + avgSmodulus += this->shearModulus(elemIdx); + } + std::size_t numDof = this->model().numGridDof(); + const auto& comm = this->simulator().vanguard().grid().comm(); + avgSmodulus = comm.sum(avgSmodulus); + Scalar numTotalDof = comm.sum(numDof); + avgSmodulus /= numTotalDof; + avgSmodulus = std::sqrt(avgSmodulus); + + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + if (eqIdx < contiRotEqIdx) { + geoMechModel_.setEqWeight(eqIdx, 1 / avgSmodulus); + } + else { + geoMechModel_.setEqWeight(eqIdx, avgSmodulus); + } + } + } + + /*! + * \brief Called by the simulator before each time integration. + */ + void beginTimeStep() override + { + // Call parent class beginTimeStep() + ParentType::beginTimeStep(); + + // Update mechanics boundary conditions. + // NOTE: Flow boundary conditions should be updated in ParentType::beginTimeStep() + if (this->nonTrivialBoundaryConditions()) { + geoMechModel_.linearizer().updateBoundaryConditionData(); + } + } + + /*! + * \brief Organize mechanics boundary conditions + * + * \param globalSpaceIdx Cell index + * \param directionId Direction id + * \returns A pair of BCMECHType and displacement vector + * + * Output from this function is used in LocalResidual::computeBoundaryTerm + * + * \note Only BCMECHTYPE = FREE and NONE implemented. FIXED will/should throw an error when computed in local + * residual! + */ + std::pair> + mechBoundaryCondition(const unsigned int globalSpaceIdx, const int directionId) + { + // Default boundary conditions if BCCON/BCPROP not defined + if (!this->nonTrivialBoundaryConditions()) { + return { BCMECHType::NONE, Dune::FieldVector{0.0, 0.0, 0.0} }; + } + + // Default for BCPROP index = 0 or no BCPROP defined at current episode + FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId); + const auto& schedule = this->simulator().vanguard().schedule(); + if (this->bcindex_(dir)[globalSpaceIdx] == 0 || schedule[this->episodeIndex()].bcprop.size() == 0) { + return { BCMECHType::NONE, Dune::FieldVector{0.0, 0.0, 0.0} }; + } + + // Get current BC + const auto& bc = schedule[this->episodeIndex()].bcprop[this->bcindex_(dir)[globalSpaceIdx]]; + if (bc.bcmechtype == BCMECHType::FREE) { + return { BCMECHType::FREE, Dune::FieldVector{0.0, 0.0, 0.0} }; + } + else { + return { bc.bcmechtype, Dune::FieldVector{0.0, 0.0, 0.0} }; + } + } + + /*! + * \brief Set mechanics source term, in particular coupling terms + * + * \param sourceTerm Source term vector + * \param globalSpaceIdx Cell index + * \param timeIdx Time index + * + * This term requires that intensive quantities are updated! + */ + void tpsaSource(Dune::FieldVector& sourceTerm, + unsigned globalSpaceIdx, + unsigned timeIdx) + { + sourceTerm = 0.0; + + // Coupling term Flow -> TPSA + const auto biot = this->biotCoeff(globalSpaceIdx); + const auto lameParam = this->lame(globalSpaceIdx); + + const auto& iq = this->model().intensiveQuantities(globalSpaceIdx, timeIdx); + const auto& fs = iq.fluidState(); + const auto pres = decay(fs.pressure(this->refPressurePhaseIdx_())); + const auto initPres = this->initialFluidState(globalSpaceIdx).pressure(this->refPressurePhaseIdx_()); + + auto sourceFromFlow = -biot / lameParam * (pres - initPres); + sourceTerm[contiSolidPresEqIdx] += sourceFromFlow; + } + + /*! + * \brief Pore volume change due to geomechanics + * + * \param elementIdx Cell index + * \param timeIdx Time index + * \returns Pore volume change (dimensionless) + * + * \note This is the coupling term to Flow + */ + Scalar rockMechPoroChange(unsigned elementIdx, unsigned timeIdx) const + { + // TODO: get timeIdx=1 solid pressure from a cached materialState (or intensiveQuantities) if/when implemented + assert (timeIdx <= historySize); + const auto solidPres = (timeIdx == 0) ? + decay( geoMechModel_.materialState(elementIdx, /*timeIdx=*/timeIdx).solidPressure()) : + geoMechModel_.solution(/*timeIdx=*/timeIdx)[elementIdx][solidPres0Idx]; + const auto biot = this->biotCoeff(elementIdx); + const auto lameParam = this->lame(elementIdx); + + return biot / lameParam * solidPres; + } + + // /// + // Public get functions + // /// + /*! + * \brief Direct access to average (half-)weight at interface between two elements + * + * \param globalElemIdxIn Inside cell index + * \param globalElemIdxOut Outside cell index + * \returns Weight average + */ + Scalar weightAverage(unsigned globalElemIdxIn, unsigned globalElemIdxOut) + { + return faceProps_.weightAverage(globalElemIdxIn, globalElemIdxOut); + } + + /*! + * \brief Direct access to normal distance at the boundary + * + * \param globalElemIdxIn Inside cell index + * \param boundaryFaceIdx Boundary (local) face index + * \returns Weight average at boundary + */ + Scalar weightAverageBoundary(unsigned globalElemIdxIn, unsigned boundaryFaceIdx) const + { + return faceProps_.weightAverageBoundary(globalElemIdxIn, boundaryFaceIdx); + } + + /*! + * \brief Direct access to product of weights at interface between two elements + * + * \param globalElemIdxIn Inside cell index + * \param globalElemIdxOut Outside cell index + * \returns Weight product + */ + Scalar weightProduct(unsigned globalElemIdxIn, unsigned globalElemIdxOut) const + { + return faceProps_.weightProduct(globalElemIdxIn, globalElemIdxOut); + } + + /*! + * \brief Direct access to normal distance between two elements + * + * \param globalElemIdxIn Inside cell index + * \param globalElemIdxOut Outside cell index + * \returns Normal distance + */ + Scalar normalDistance(unsigned globalElemIdxIn, unsigned globalElemIdxOut) const + { + return faceProps_.normalDistance(globalElemIdxIn, globalElemIdxOut); + } + + /*! + * \brief Direct access to normal distance at the boundary + * + * \param globalElemIdxIn Inside cell index + * \param boundaryFaceIdx Boundary (local) face index + * \returns Normal distance to boundary + */ + Scalar normalDistanceBoundary(unsigned globalElemIdxIn, unsigned boundaryFaceIdx) const + { + return faceProps_.normalDistanceBoundary(globalElemIdxIn, boundaryFaceIdx); + } + + /*! + * \brief Direct access to face normal between two elements + * + * \param globalElemIdxIn Inside cell index + * \param globalElemIdxOut Outside cell index + * \returns Face normal + */ + DimVector cellFaceNormal(unsigned globalElemIdxIn, unsigned globalElemIdxOut) + { + return faceProps_.cellFaceNormal(globalElemIdxIn, globalElemIdxOut); + } + + /*! + * \brief Direct access to face normal at the boundary + * + * \param globalElemIdxIn Inside cell index + * \param boundaryFaceIdx Boundary (local) face index + * \returns Face normal on boundary + */ + const DimVector& cellFaceNormalBoundary(unsigned globalElemIdxIn, unsigned boundaryFaceIdx) const + { + return faceProps_.cellFaceNormalBoundary(globalElemIdxIn, boundaryFaceIdx); + } + + /*! + * \brief Direct access to shear modulus in an element + * + * \param globalElemIdx Cell index + * \returns Shear modulus + */ + Scalar shearModulus(unsigned globalElemIdx) const + { + return faceProps_.shearModulus(globalElemIdx); + } + + /*! + * \brief Flow-TPSA lagged coupling scheme activated? + * + * \returns Bool indicating lagged coupling + */ + bool laggedScheme() const + { + const auto& mechSolver = this->simulator().vanguard().eclState().runspec().mechSolver(); + return mechSolver.laggedScheme(); + } + + /*! + * \brief Flow-TPSA fixed-stress coupling scheme activated? + * + * \returns Bool indicating fixed-stress coupling + */ + bool fixedStressScheme() const + { + const auto& mechSolver = this->simulator().vanguard().eclState().runspec().mechSolver(); + return mechSolver.fixedStressScheme(); + } + + /*! + * \brief Get TPSA model + * + * \returns Reference to geomechanics model + */ + const GeomechModel& geoMechModel() const + { + return geoMechModel_; + } + + /*! + * \brief Get TPSA model + * + * \returns Reference to geomechanics model + */ + GeomechModel& geoMechModel() + { + return geoMechModel_; + } + + /*! + * \brief Get fixed-stress iteration parameters + * + * \returns Pair with min/max fixed-stress iterations + */ + std::pair fixedStressParameters() const + { + const auto& mechSolver = this->simulator().vanguard().eclState().runspec().mechSolver(); + return std::make_pair(mechSolver.fixedStressMinIter(), mechSolver.fixedStressMaxIter()); + } + +protected: + // /// + // Protected functions + // /// + /*! + * \brief Read initial conditions and generate material state for TPSA model + */ + void readInitalConditionsTPSA_() + { + // /// + // OBS: No equilibration keywords (e.g., STREQUIL) implemented yet! + // /// + + // Set all initial material state variables to zero + std::size_t numDof = this->model().numGridDof(); + initialMaterialState_.resize(numDof); + for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + auto& dofMaterialState = initialMaterialState_[dofIdx]; + for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) { + dofMaterialState.setDisplacement(dirIdx, 0.0); + dofMaterialState.setRotation(dirIdx, 0.0); + } + dofMaterialState.setSolidPressure(0.0); + } + } + +private: + FaceProperties faceProps_; + GeomechModel geoMechModel_; + + std::vector biotcoeff_; + std::vector initialMaterialState_; +}; // class FlowProblemTPSA + +} // namespace Opm + +#endif diff --git a/opm/simulators/flow/MainDispatchDynamic.cpp b/opm/simulators/flow/MainDispatchDynamic.cpp index 2e8a459726..e1ee8f0e73 100644 --- a/opm/simulators/flow/MainDispatchDynamic.cpp +++ b/opm/simulators/flow/MainDispatchDynamic.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -43,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -54,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -178,6 +181,8 @@ int Opm::Main::runTwoPhase(const Phases& phases) const bool disgasw = eclipseState_->getSimulationConfig().hasDISGASW(); const bool vapwat = eclipseState_->getSimulationConfig().hasVAPWAT(); + const auto& rspec = this->eclipseState_->runspec(); + // oil-gas if (phases.active(Phase::OIL) && phases.active(Phase::GAS)) { if (diffusive) { @@ -208,6 +213,8 @@ int Opm::Main::runTwoPhase(const Phases& phases) return flowGasWaterDissolutionDiffuseMain(argc_, argv_, outputCout_, outputFiles_); + } else if (rspec.mech() && rspec.mechSolver().tpsa()) { + return flowGasWaterDissolutionTpsaMain(argc_, argv_, outputCout_, outputFiles_); } return flowGasWaterDissolutionMain(argc_, argv_, outputCout_, outputFiles_); @@ -286,6 +293,8 @@ int Opm::Main::runFoam() int Opm::Main::runWaterOnly(const Phases& phases) { + const auto& rspec = this->eclipseState_->runspec(); + if (!phases.active(Phase::WATER) || phases.size() != 1) { if (outputCout_) { std::cerr << "No valid configuration is found for " @@ -296,6 +305,10 @@ int Opm::Main::runWaterOnly(const Phases& phases) return EXIT_FAILURE; } + if (rspec.mech() && rspec.mechSolver().tpsa()) { + return flowWaterOnlyTpsaMain(argc_, argv_, outputCout_, outputFiles_); + } + return flowWaterOnlyMain(argc_, argv_, outputCout_, outputFiles_); } @@ -438,6 +451,12 @@ int Opm::Main::runBlackOil() return flowBlackoilMain(argc_, argv_, outputCout_, outputFiles_); } if (this->eclipseState_->runspec().hysterPar().active()) { + const auto& rspec = this->eclipseState_->runspec(); + if (rspec.mech() && rspec.mechSolver().tpsa()) { + // Blackoil + TPSA geomechanics + return flowBlackoilTpsaMain(argc_, argv_, outputCout_, outputFiles_); + } + return flowBlackoilTpfaMain(argc_, argv_, outputCout_, outputFiles_); } else { // Use variant without hysteresis support to save memory. diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 86c029c2d4..233dda1fad 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -38,6 +38,8 @@ #include +#include + #include #include #include @@ -47,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -159,6 +162,9 @@ public: SolverParameters::registerParameters(); TimeStepper::registerParameters(); detail::registerSimulatorParameters(); + + TpsaNewtonMethodParams::registerParameters(); + TpsaLinearSolverParameters::registerParameters(); } /// Run the simulation. diff --git a/opm/simulators/flow/TTagFlowProblemTPSA.hpp b/opm/simulators/flow/TTagFlowProblemTPSA.hpp new file mode 100644 index 0000000000..45eed41f99 --- /dev/null +++ b/opm/simulators/flow/TTagFlowProblemTPSA.hpp @@ -0,0 +1,150 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TTAG_FLOW_PROBLEM_TPSA_HPP +#define TTAG_FLOW_PROBLEM_TPSA_HPP + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +namespace Opm::Properties { + +namespace TTag { + +struct FlowProblemTpsa {}; + +} // Opm::Properties::TTag + +// TPSA indices for primary variables and equations +template +struct IndicesTPSA +{ + using type = ElasticityIndices; +}; + +// Number of TPSA equations +template +struct NumEqTPSA +{ static constexpr int value = GetPropType::numEq; }; + +// TPSA linearizer +template +struct LinearizerTPSA +{ using type = TpsaLinearizer; }; + +// Set the function evaluation w.r.t. the TPSA primary variables +template +struct EvaluationTPSA +{ +private: + static constexpr unsigned numEq = getPropValue(); + + using Scalar = GetPropType; + +public: + using type = DenseAd::Evaluation; +}; + +// TPSA equation vector +template +struct EqVectorTPSA +{ + using type = Dune::FieldVector, + getPropValue()>; +}; + +// Global TPSA equation vector +template +struct GlobalEqVectorTPSA +{ using type = Dune::BlockVector>; }; + +// TPSA Newton method +template +struct NewtonMethodTPSA +{ using type = TpsaNewtonMethod; }; + +// TPSA primary variables +template +struct PrimaryVariablesTPSA +{ using type = ElasticityPrimaryVariables; }; + +// TPSA solution vector +template +struct SolutionVectorTPSA +{ using type = Dune::BlockVector>; }; + +// TPSA number of historic solutions to save +template +struct SolutionHistorySizeTPSA +{ static constexpr int value = 2; }; + +// TPSA model +template +struct ModelTPSA +{ using type = TpsaModel; }; + +// TPSA local residual +template +struct LocalResidualTPSA +{ using type = ElasticityLocalResidual; }; + +// TPSA sparse matrix adapter for Jacobian +template +struct SparseMatrixAdapterTPSA +{ +private: + using Scalar = GetPropType; + enum { numEq = getPropValue() }; + using Block = MatrixBlock; + +public: + using type = typename Linear::IstlSparseMatrixAdapter; + +}; + +// Set linear solver backend +template +struct LinearSolverBackendTPSA +{ using type = ISTLSolverTPSA; }; + +} // namespace Opm::Properties + +#endif diff --git a/opm/simulators/linalg/AbstractISTLSolver.hpp b/opm/simulators/linalg/AbstractISTLSolver.hpp index ba6dbee98a..caf1ce29a6 100644 --- a/opm/simulators/linalg/AbstractISTLSolver.hpp +++ b/opm/simulators/linalg/AbstractISTLSolver.hpp @@ -40,7 +40,7 @@ namespace Opm * directly. This would simplify the interface and reduce the number of * required method calls before solving the system. */ -template +template class AbstractISTLSolver { public: @@ -50,8 +50,6 @@ public: using CommunicationType = Dune::Communication; #endif - using SparseMatrixAdapter = GetPropType; - using Vector = GetPropType; using Matrix = typename SparseMatrixAdapter::IstlMatrix; virtual ~AbstractISTLSolver() = default; diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index e55762e8d3..8ecb0e4be1 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -145,7 +145,8 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, /// as a block-structured matrix (one block for all cell variables) for a fixed /// number of cell variables np . template - class ISTLSolver : public AbstractISTLSolver + class ISTLSolver : public AbstractISTLSolver, + GetPropType> { protected: using GridView = GetPropType; @@ -464,7 +465,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, bool checkConvergence(const Dune::InverseOperatorResult& result) const { - return AbstractISTLSolver::checkConvergence(result, parameters_[activeSolverNum_]); + return AbstractISTLSolver::checkConvergence(result, parameters_[activeSolverNum_]); } bool isParallel() const { diff --git a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp index ee9844262b..9be9e5796f 100644 --- a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp +++ b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp @@ -40,7 +40,8 @@ namespace Opm * \brief ISTLSolverRuntimeOptionProxy selects the appropriate ISTLSolver runtime based on CLI options */ template -class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver +class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver, + GetPropType> { public: using SparseMatrixAdapter = GetPropType; @@ -143,7 +144,7 @@ public: } private: - std::unique_ptr> istlSolver_; + std::unique_ptr> istlSolver_; template diff --git a/opm/simulators/linalg/ISTLSolverTPSA.hpp b/opm/simulators/linalg/ISTLSolverTPSA.hpp new file mode 100644 index 0000000000..e3e86094a0 --- /dev/null +++ b/opm/simulators/linalg/ISTLSolverTPSA.hpp @@ -0,0 +1,415 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef ISTL_SOLVER_TPSA_HPP +#define ISTL_SOLVER_TPSA_HPP + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +namespace Opm { + +/*! +* \brief Class for setting up ISTL linear solvers for TPSA +* +* Implements AbstractISTLSolver interface for TPSA linear solvers. +*/ +template +class ISTLSolverTPSA : public AbstractISTLSolver, + GetPropType> +{ + using ElementMapper = GetPropType; + using Simulator = GetPropType; + using SparseMatrixAdapter = GetPropType; + using Vector = GetPropType; + + using Matrix = typename SparseMatrixAdapter::IstlMatrix; + + +#if HAVE_MPI + using CommunicationType = Dune::OwnerOverlapCopyCommunication; +#else + using CommunicationType = Dune::Communication; +#endif + +public: + /*! + * \brief Constructor + * + * \param simulator Simulator object + */ + ISTLSolverTPSA(const Simulator& simulator) + : simulator_(simulator) + , solveCount_(0) + , iterations_(0) + , matrix_(nullptr) + , rhs_(nullptr) + { + // Init parameters + parameters_.init(); + + // Initialize linear solver + initialize(); + } + + /*! + * \brief Register runtime/default parameters for linear solver + */ + static void registerParameters() + { + TpsaLinearSolverParameters::registerParameters(); + } + + /*! + * \brief Setup linear solver object based on runtime/default parameters + */ + void initialize() + { + // Setup property tree for FlexibleSolver + prm_ = setupPropertyTree(parameters_, false, false, /*tpsaSetup=*/true); + + // Reset comm_ pointer +#if HAVE_MPI + comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); +#endif + + // Extract and copy parallel grid information + extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); +#if HAVE_MPI + if (isParallel()) { + const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0); + detail::copyParValues(parallelInformation_, size, *comm_); + } +#endif + + // Get info on overlapping rows + ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout()); + std::vector dummyInteriorRows; + detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, dummyInteriorRows); + + // Set number of interior cells in FlexibleSolverInfo + flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true); + + // Print parameters to PRT/DBG logs. + detail::printLinearSolverParameters(parameters_, prm_, simulator_.gridView().comm()); + } + + /*! + * \brief Prepare matix and rhs vector for linear solve + * + * \param M System matrix + * \param b Right-hand side vector + */ + void initPrepare(const Matrix& M, Vector& b) + { + // Update matrix entries if it has not been initialized yet + const bool firstcall = (matrix_ == nullptr); + if (firstcall) { + // Model will not change the matrix object. Hence simply store a pointer to the original one with a deleter + // that does nothing. + // OBS: We need to be able to scale the linear system, hence const_cast + matrix_ = const_cast(&M); + } + else { + // Pointers should not change; throw if the case + if (&M != matrix_) { + OPM_THROW(std::logic_error, "TPSA: Matrix objects are expected to be reused when reassembling!"); + } + + } + + // Set right-hand side vector + rhs_ = &b; + + // Zero out the overlapping cells in matrix (not in ilu0 case) + std::string type = prm_.template get("preconditioner.type", "paroverilu0"); + std::transform(type.begin(), type.end(), type.begin(), ::tolower); + if (isParallel() && type != "paroverilu0") { + detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_); + } + } + + /*! + * \copydoc AbstractISTLSolver::prepare + */ + void prepare(const Matrix& M, Vector& b) override + { + try { + initPrepare(M, b); + + prepareFlexibleSolver(); + } + OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR + ("TPSA: Failure likely due to a faulty linear solver JSON specification. " + "Check for errors related to missing nodes."); + } + + /*! + * \copydoc AbstractISTLSolver::prepare + */ + void prepare(const SparseMatrixAdapter& M, Vector& b) override + { + prepare(M.istlMatrix(), b); + } + + /*! + * \brief Prepare linear solver + * + * Create linear solver using FlexibleSolverInfo struct from ISTLSolver.hpp, or update preconditioner if solver has + * been created + */ + void prepareFlexibleSolver() + { + // Create solver or just update preconditioner + if (!flexibleSolver_.solver_) { + // Dummy weights calculator + // TODO: TPSA specific weights calculation for AMG preconditioner + std::function weightCalculator; + + // Create FlexibleSolver + flexibleSolver_.create(getMatrix(), + isParallel(), + prm_, + /*pressureIndex=*/0, + weightCalculator, + /*forceSerial_=*/false, + comm_.get()); + } + else { + // Update preconditioner + flexibleSolver_.pre_->update(); + } + } + + /*! + * \copydoc AbstractISTLSolver::solve + */ + bool solve(Vector& x) override + { + // Increase solver count + ++solveCount_; + + // Write system matrix if verbosity level is high + const int verbosity = prm_.get("verbosity", 0); + if (verbosity > 10) { + // simulator_ is only used to get names + Helper::writeSystem(simulator_, + getMatrix(), + *rhs_, + comm_.get()); + } + + // Solve linear system + Dune::InverseOperatorResult result; + assert(flexibleSolver_.solver_); + flexibleSolver_.solver_->apply(x, *rhs_, result); + + // Store no. linear iterations + iterations_ = result.iterations; + + // Return result for convergence check (boolean) + return checkConvergence(result); + } + + /*! + * \brief Reset number of solver calls to zero + */ + void resetSolveCount() + { + solveCount_ = 0; + } + + // //// + // Unused AbstractISTLSolver functions + // //// + /*! + * \copydoc AbstractISTLSolver::eraseMatrix + */ + void eraseMatrix() override + { } + + /*! + * \copydoc AbstractISTLSolver::setActiveSolver + */ + void setActiveSolver(int /*num*/) override + { } + + /*! + * \copydoc AbstractISTLSolver::numAvailableSolvers + */ + int numAvailableSolvers() const override + { + return 1; + } + + /*! + * \copydoc AbstractISTLSolver::setResidual + */ + void setResidual(Vector& /*b*/) override + { } + + /*! + * \copydoc AbstractISTLSolver::setMatrix + */ + void setMatrix(const SparseMatrixAdapter& /*M*/) override + { } + + // /// + // Public get functions + // /// + /*! + * \copydoc AbstractISTLSolver::getResidual + */ + void getResidual(Vector& b) const override + { + b = *rhs_; + } + + /*! + * \copydoc AbstractISTLSolver::getSolveCount + */ + int getSolveCount() const override + { + return solveCount_; + } + + /*! + * \copydoc AbstractISTLSolver::iterations + */ + int iterations() const override + { + return iterations_; + } + + /*! + * \copydoc AbstractISTLSolver::comm + */ + const CommunicationType* comm() const override + { + return comm_.get(); + } + +protected: + /*! + * \brief Check for parallel session + * + * \returns Bool indicating if this is a parallel session + * + * \warning comm_ must be set before using this function + */ + bool isParallel() const + { +#if HAVE_MPI + return comm_->communicator().size() > 1; +#else + return false; +#endif + } + + /*! + * \brief Check for linear solver convergence + * + * \param result Linear solver result container + * \returns Bool indicating convergence + */ + bool checkConvergence(const Dune::InverseOperatorResult& result) const + { + return AbstractISTLSolver::checkConvergence(result, parameters_); + } + + // /// + // Protected get functions + // /// + /*! + * \brief Get reference to system matrix object + * + * \returns Reference to matrix + */ + Matrix& getMatrix() + { + if (!matrix_) { + OPM_THROW(std::runtime_error, "TPSA: System matrix \"M\" not defined!"); + } + return *matrix_; + } + + /*! + * \brief Get reference to system matrix object + * + * \returns Reference to matrix + */ + const Matrix& getMatrix() const + { + if (!matrix_) { + OPM_THROW(std::runtime_error, "TPSA: System matrix \"M\" not defined!"); + } + return *matrix_; + } + + const Simulator& simulator_; + std::any parallelInformation_; + int solveCount_; + int iterations_; + + detail::FlexibleSolverInfo flexibleSolver_; + TpsaLinearSolverParameters parameters_; + PropertyTree prm_; + + Matrix* matrix_; + Vector* rhs_; + + std::shared_ptr comm_; + std::vector overlapRows_; +}; + +} // namespace Opm + +#endif diff --git a/opm/simulators/linalg/TPSALinearSolverParameters.cpp b/opm/simulators/linalg/TPSALinearSolverParameters.cpp new file mode 100644 index 0000000000..5734924608 --- /dev/null +++ b/opm/simulators/linalg/TPSALinearSolverParameters.cpp @@ -0,0 +1,108 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include + +#include + + +namespace Opm { + +/*! +* \brief Internalize runtime parameters +* +* Overloading FlowLinearSolverParameters::init() to read TPSA specific runtime/default parameters +*/ +void TpsaLinearSolverParameters::init() +{ + // Runtime parameters + linear_solver_reduction_ = Parameters::Get(); + relaxed_linear_solver_reduction_ = Parameters::Get(); + linear_solver_maxiter_ = Parameters::Get(); + linear_solver_restart_ = Parameters::Get(); + linear_solver_verbosity_ = Parameters::Get(); + ilu_relaxation_ = Parameters::Get(); + ilu_fillin_level_ = Parameters::Get(); + newton_use_gmres_ = Parameters::Get(); + ignoreConvergenceFailure_ = Parameters::Get(); + linsolver_ = Parameters::Get(); + linear_solver_print_json_definition_ = Parameters::Get(); + + // Hardcode use of CPU linear solvers (?) + linear_solver_accelerator_ = Parameters::LinearSolverAcceleratorType::CPU; +} + +/*! +* \brief Register TPSA linear solver runtime parameters +*/ +void TpsaLinearSolverParameters::registerParameters() +{ + Parameters::Register + ("Minimum residual reduction in TPSA linear solver for convergenc"); + Parameters::Register + ("A relaxed version of --tpsa-linear-solver-reduction (use with care!)"); + Parameters::Register + ("Maximum TPSA linear iterations"); + Parameters::Register + ("Number of iterations before restarting GMRES if --tpsa-use-gmres=true"); + Parameters::Register + ("Level of verbosity in TPSA linear solver: 0 = off, 2 = all"); + Parameters::Register + ("Relaxation factor for TPSA linear solver ILU preconditioner"); + Parameters::Register + ("Fill-in level of TPSA linear solver ILU preconditioner"); + Parameters::Register + ("Use GMRES linear solver. If false, BiCGStab is used."); + Parameters::Register + ("Continue simulation even if TPSA linear solver did not converge"); + Parameters::Register + ("Configuration for linear solver. Valid preset options are: ilu0, dilu, amg or umfpack. " + "Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, " + "ending with '.json.'"); + Parameters::Register + ("Print JSON formatted configuration of the TPSA linear solver. Can be used to make configuration JSON file " + "for --tpsa-linear-solver"); +} + +/*! +* \brief Reset TPSA linear solver parameters +* +* \warning Should be the same as the corresponding Parameters defaults! +*/ +void TpsaLinearSolverParameters::reset() +{ + linear_solver_reduction_ = 1e-3; + relaxed_linear_solver_reduction_ = 1e-3; + linear_solver_maxiter_ = 200; + linear_solver_restart_ = 40; + linear_solver_verbosity_ = 0; + ilu_relaxation_ = 9; + ilu_fillin_level_ = 0; + newton_use_gmres_ = false; + ignoreConvergenceFailure_ = false; + linsolver_ = "ilu0"; + linear_solver_print_json_definition_ = false; +} + +} // namespace Opm diff --git a/opm/simulators/linalg/TPSALinearSolverParameters.hpp b/opm/simulators/linalg/TPSALinearSolverParameters.hpp new file mode 100644 index 0000000000..fb43adea32 --- /dev/null +++ b/opm/simulators/linalg/TPSALinearSolverParameters.hpp @@ -0,0 +1,64 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef TPSA_LINEAR_SOLVER_PARAMETERS_HPP +#define TPSA_LINEAR_SOLVER_PARAMETERS_HPP + +#include + +#include + + +// Default runtime parameters +namespace Opm::Parameters { + +struct TpsaLinearSolverReduction { static constexpr double value = 1e-3; }; +struct TpsaRelaxedLinearSolverReduction { static constexpr double value = 1e-3; }; +struct TpsaLinearSolverMaxIter { static constexpr int value = 200; }; +struct TpsaLinearSolverRestart { static constexpr int value = 40; }; +struct TpsaLinearSolverVerbosity { static constexpr int value = 0; }; +struct TpsaIluRelaxation { static constexpr double value = 0.9; }; +struct TpsaIluFillinLevel { static constexpr int value = 0; }; +struct TpsaUseGmres { static constexpr bool value = false; }; +struct TpsaLinearSolverIgnoreConvergenceFailure { static constexpr bool value = false; }; +struct TpsaLinearSolver { static constexpr auto value = "ilu0"; }; +struct TpsaLinearSolverPrintJsonDefinition { static constexpr auto value = false; }; + +} // namespace Opm::Parameters + +namespace Opm { + +/*! +* \brief Parametern for linear solver and preconditioner +*/ +struct TpsaLinearSolverParameters : public FlowLinearSolverParameters +{ + void init(); + static void registerParameters(); + void reset(); +}; + +} // namespace Opm + +#endif diff --git a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp index 68c89c86e3..d4f4d012cb 100644 --- a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp +++ b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp @@ -57,7 +57,8 @@ namespace Opm::gpuistl * matrices and vectors internally for computations. */ template -class ISTLSolverGPUISTL : public AbstractISTLSolver +class ISTLSolverGPUISTL : public AbstractISTLSolver, + GetPropType> { public: using SparseMatrixAdapter = GetPropType; @@ -332,7 +333,7 @@ public: private: bool checkConvergence(const Dune::InverseOperatorResult& result) const { - return AbstractISTLSolver::checkConvergence(result, m_parameters); + return AbstractISTLSolver::checkConvergence(result, m_parameters); } // Weights to make approximate pressure equations. diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 3f7b6919ec..6833ca20aa 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -27,6 +27,7 @@ #include #include +#include namespace Opm { @@ -180,7 +181,8 @@ namespace PropertyTree setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters to potentially override. bool linearSolverMaxIterSet, - bool linearSolverReductionSet) + bool linearSolverReductionSet, + bool tpsaSetup) { std::string conf = p.linsolver_; @@ -201,28 +203,30 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters std::ranges::transform(conf, conf.begin(), ::tolower); // Use CPR configuration. - if ((conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes") || (conf == "cpr_trueimpesanalytic")) { - if (!linearSolverMaxIterSet) { - // Use our own default unless it was explicitly overridden by user. - p.linear_solver_maxiter_ = 20; + if (!tpsaSetup) { + if ((conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes") || (conf == "cpr_trueimpesanalytic")) { + if (!linearSolverMaxIterSet) { + // Use our own default unless it was explicitly overridden by user. + p.linear_solver_maxiter_ = 20; + } + if (!linearSolverReductionSet) { + // Use our own default unless it was explicitly overridden by user. + p.linear_solver_reduction_ = 0.005; + } + return setupCPR(conf, p); } - if (!linearSolverReductionSet) { - // Use our own default unless it was explicitly overridden by user. - p.linear_solver_reduction_ = 0.005; - } - return setupCPR(conf, p); - } - if ((conf == "cpr") || (conf == "cprw")) { - if (!linearSolverMaxIterSet) { - // Use our own default unless it was explicitly overridden by user. - p.linear_solver_maxiter_ = 20; + if ((conf == "cpr") || (conf == "cprw")) { + if (!linearSolverMaxIterSet) { + // Use our own default unless it was explicitly overridden by user. + p.linear_solver_maxiter_ = 20; + } + if (!linearSolverReductionSet) { + // Use our own default unless it was explicitly overridden by user. + p.linear_solver_reduction_ = 0.005; + } + return setupCPRW(conf, p); } - if (!linearSolverReductionSet) { - // Use our own default unless it was explicitly overridden by user. - p.linear_solver_reduction_ = 0.005; - } - return setupCPRW(conf, p); } if (conf == "amg") { @@ -261,9 +265,16 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters } // No valid configuration option found. - OPM_THROW(std::invalid_argument, - conf + " is not a valid setting for --linear-solver-configuration." - " Please use ilu0, dilu, cpr, cprw, cpr_trueimpes, cpr_quasiimpes, cpr_trueimpesanalytic or isai"); + if (tpsaSetup) { + OPM_THROW(std::invalid_argument, + fmt::format("No valid settings found for --tpsa-linear-solver={}! " + "Valid preset options are: ilu0, dilu, amg, or umfpack.", conf)); + } + else { + OPM_THROW(std::invalid_argument, + conf + " is not a valid setting for --linear-solver-configuration." + " Please use ilu0, dilu, cpr, cprw, cpr_trueimpes, cpr_quasiimpes, cpr_trueimpesanalytic or isai"); + } } std::string getSolverString(const FlowLinearSolverParameters& p) diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp index 2cb1cc20f2..a424a57325 100644 --- a/opm/simulators/linalg/setupPropertyTree.hpp +++ b/opm/simulators/linalg/setupPropertyTree.hpp @@ -31,7 +31,8 @@ struct FlowLinearSolverParameters; PropertyTree setupPropertyTree(FlowLinearSolverParameters p, bool linearSolverMaxIterSet, - bool linearSolverReductionSet); + bool linearSolverReductionSet, + bool tpsaSetup = false); PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters& p); PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p); diff --git a/regressionTests.cmake b/regressionTests.cmake index f8d30b2e3f..c395e334ed 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -639,6 +639,21 @@ add_multiple_tests( TEST_ARGS --tolerance-cnv-relaxed=0.01 ) +set(_tpsa_cases + TPSA_LAGGED + TPSA_FIXEDSTRESS +) + +add_multiple_tests( + _tpsa_cases + "" + SIMULATOR flow + ABS_TOL ${abs_tol} + REL_TOL ${rel_tol} + DIR tpsa + TEST_ARGS --enable-opm-rst-file=1 +) + add_test_compareECLFiles(CASENAME ppcwmax FILENAME PPCWMAX-01 SIMULATOR flow diff --git a/tests/test_tpsa_face_properties.cpp b/tests/test_tpsa_face_properties.cpp new file mode 100644 index 0000000000..86056602d7 --- /dev/null +++ b/tests/test_tpsa_face_properties.cpp @@ -0,0 +1,192 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025 NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include + +#define BOOST_TEST_MODULE TestFacePropertiesTPSA +#define BOOST_TEST_NO_MAIN + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#if HAVE_MPI +struct MPIError +{ + MPIError(std::string s, int e) : errorstring(std::move(s)), errorcode(e){} + std::string errorstring; + int errorcode; +}; + +void MPI_err_handler(MPI_Comm*, int* err_code, ...) +{ + std::vector err_string(MPI_MAX_ERROR_STRING); + int err_length; + MPI_Error_string(*err_code, err_string.data(), &err_length); + std::string s(err_string.data(), err_length); + std::cerr << "An MPI Error ocurred:" << std::endl << s << std::endl; + throw MPIError(s, *err_code); +} +#endif + +bool init_unit_test_func() +{ + return true; +} + +static Opm::Deck createDeck() +{ + // Deck to test + std::string deck_string = R"( + RUNSPEC + DIMENS + 2 2 2 / + + GRID + DX + 8*50.0 / + DY + 8*50.0 / + DZ + 4*30.0 4*50.0 / + TOPS + 3*1000.0 1*1020 / + + PORO + 8*0.3 / + SMODULUS + 8*1.0 / + + END + )"; + Opm::Parser parser; + return parser.parseString(deck_string); +} + +BOOST_AUTO_TEST_CASE(SimpleGridWithNNC) +{ + // using declarations + using Grid = Dune::CpGrid; + using GridView = Grid::LeafGridView; + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + using CartesianIndexMapper = Dune::CartesianIndexMapper; + using FacePropertiesTPSA = Opm::FacePropertiesTPSA; + using DimVector = Dune::FieldVector; + + // Test simple 2x2x2 grid with one column shifted downwards + Opm::Deck deck = createDeck(); + Opm::EclipseState eclState(deck); + + Grid grid; + grid.processEclipseFormat(&eclState.getInputGrid(), &eclState, false, false, false); + const auto& gridView = grid.leafGridView(); + + CartesianIndexMapper cartMapper = Dune::CartesianIndexMapper(grid); + auto centroids = [&eclState, &cartMapper](int index) + { return eclState.getInputGrid().getCellCenter(cartMapper.cartesianIndex(index)); }; + + // Init. FacePropertiesTPSA and calculate all properties + FacePropertiesTPSA faceProps(eclState, + gridView, + cartMapper, + grid, + centroids); + faceProps.update(); + + // Check face properties + // Natural neighbor checks between cells 0 and 4 + const unsigned elem1 = 0; + const unsigned elem2 = 4; + const double normDist = 40.0; + const double weightAvg_04 = 0.375; + const double weightAvg_40 = 0.625; // = 1 - weightAvg_04 + const double weightProd = 3.75e-16; + const DimVector normal_04 = {0.0, 0.0, 1.0}; + const DimVector normal_40 = {0.0, 0.0, -1.0}; // = -1 * normal_04 + BOOST_CHECK_CLOSE(faceProps.weightAverage(elem1, elem2), weightAvg_04, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightAverage(elem2, elem1), weightAvg_40, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightProduct(elem1, elem2), weightProd, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.normalDistance(elem1, elem2), normDist, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.normalDistance(elem2, elem1), normDist, 1.0e-8); + BOOST_CHECK_EQUAL(faceProps.cellFaceNormal(elem1, elem2), normal_04); + BOOST_CHECK_EQUAL(faceProps.cellFaceNormal(elem2, elem1), normal_40); + + // NNC checks between cell 3 and 5 + const unsigned elemNNC1 = 3; + const unsigned elemNNC2 = 5; + const double normDistNNC = 50.0; + const double weightAvgNNC = 0.5; + const double weightProdNNC = 6.25e-16; + const DimVector normalNNC_35 = {0.0, -1.0, 0.0}; + const DimVector normalNNC_53 = {0.0, 1.0, 0.0}; + BOOST_CHECK_CLOSE(faceProps.weightAverage(elemNNC1, elemNNC2), weightAvgNNC, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightAverage(elemNNC2, elemNNC1), weightAvgNNC, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightProduct(elemNNC1, elemNNC2), weightProdNNC, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightProduct(elemNNC2, elemNNC1), weightProdNNC, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.normalDistance(elemNNC1, elemNNC2), normDistNNC, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.normalDistance(elemNNC2, elemNNC1), normDistNNC, 1.0e-8); + BOOST_CHECK_EQUAL(faceProps.cellFaceNormal(elemNNC1, elemNNC2), normalNNC_35); + BOOST_CHECK_EQUAL(faceProps.cellFaceNormal(elemNNC2, elemNNC1), normalNNC_53); + + // Boundary interfaces in cell 7 + const unsigned elemBnd = 7; + const double normDistBnd = 25.0; + const double weightAvgBnd = 1.0; + const double weightProdBnd = 0.0; + const std::vector normalBnd = { {-1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, // x-dir + {0.0, -1.0, 0.0}, {0.0, 1.0, 0.0}, // y-dir + {0.0, 0.0, 1.0} }; // z-dir + for (std::size_t bndIdx = 0; bndIdx < 5; ++bndIdx) { + BOOST_CHECK_CLOSE(faceProps.weightAverageBoundary(elemBnd, bndIdx), weightAvgBnd, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.weightProductBoundary(elemBnd, bndIdx), weightProdBnd, 1.0e-8); + BOOST_CHECK_CLOSE(faceProps.normalDistanceBoundary(elemBnd, bndIdx), normDistBnd, 1.0e-8); + BOOST_CHECK_EQUAL(faceProps.cellFaceNormalBoundary(elemBnd, bndIdx), normalBnd[bndIdx]); + } +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); +#if HAVE_MPI + // register a throwing error handler to allow for + // debugging with "catch throw" in gdb + MPI_Errhandler handler; + MPI_Comm_create_errhandler(MPI_err_handler, &handler); + MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler); +#endif + return boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv); +} diff --git a/tests/test_tpsa_localresidual.cpp b/tests/test_tpsa_localresidual.cpp new file mode 100644 index 0000000000..0a47367fa1 --- /dev/null +++ b/tests/test_tpsa_localresidual.cpp @@ -0,0 +1,284 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include + +#define BOOST_TEST_MODULE TpsaLocalResidualTests + +#include + +#include + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +#include +#include + +namespace Opm::Properties { + namespace TTag { + struct TpsaTestTypeTag { + using InheritsFrom = std::tuple; + }; + } + + template + struct Problem + { using type = FlowProblemTPSA; }; + + template + struct EnableMech + { static constexpr bool value = true; }; + + template + struct Indices + { + private: + using BaseTypeTag = TTag::FlowProblem; + using FluidSystem = GetPropType; + static constexpr EnergyModules energyModuleType = getPropValue(); + static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal; + + public: + using type = BlackOilOnePhaseIndices(), + getPropValue(), + getPropValue(), + numEnergyVars, + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*enabledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()>; + }; +} + +namespace { + +template +std::unique_ptr> +initSimulator(const char *filename) +{ + using Simulator = Opm::GetPropType; + const auto filenameArg = std::string {"--ecl-deck-file-name="} + filename; + const char* argv[] = { + "test_tpsa", + filenameArg.c_str() + }; + + Opm::Parameters::reset(); + Opm::registerAllParameters_(false); + Opm::BlackoilModelParameters::registerParameters(); + Opm::registerEclTimeSteppingParameters(); + Opm::Parameters::Register("Do *NOT* use!"); + Opm::Parameters::SetDefault(1); + Opm::Parameters::endRegistration(); + Opm::setupParameters_(/*argc=*/sizeof(argv) / sizeof(argv[0]), + argv, + /*registerParams=*/false, + /*allowUnused*/false, + /*handleHelp*/true, + /*myRank*/0); + + Opm::FlowGenericVanguard::readDeck(filename); + + return std::make_unique(); +} + +struct TpsaLocalResidualFixture +{ + // Constructor + TpsaLocalResidualFixture() + { + int argc = boost::unit_test::framework::master_test_suite().argc; + char** argv = boost::unit_test::framework::master_test_suite().argv; + #if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); + #else + Dune::MPIHelper::instance(argc, argv); + #endif + Opm::FlowGenericVanguard::setCommunication(std::make_unique()); + } +}; + +struct BoundaryConditionData +{ + Opm::BCMECHType type; + std::vector displacement; + unsigned boundaryFaceIndex; + double faceArea; +}; + +} // namespace Anonymous + +BOOST_GLOBAL_FIXTURE(TpsaLocalResidualFixture); + +BOOST_AUTO_TEST_CASE(TestElasticityResidual) { + using TypeTag = Opm::Properties::TTag::TpsaTestTypeTag; + using LocalResidual = Opm::ElasticityLocalResidual; + using Evaluation = Opm::GetPropType; + using MaterialState = Opm::MaterialStateTPSA; + using FluidSystem = Opm::GetPropType; + using FVector = Dune::FieldVector; + + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + + // Setup + auto simulator = initSimulator("tpsa_ex.data"); + simulator->model().applyInitialSolution(); + // simulator->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + auto& problem = simulator->problem(); + MaterialState materialStateIn; + MaterialState materialStateEx; + unsigned globI = 0; + unsigned globJ = 1; + + // Set TPSA quantities + materialStateIn.setDisplacement(0, 1.0); + materialStateIn.setDisplacement(1, 1.0); + materialStateIn.setDisplacement(2, 1.0); + materialStateIn.setRotation(0, 1.0); + materialStateIn.setRotation(1, 1.0); + materialStateIn.setRotation(2, 1.0); + materialStateIn.setSolidPressure(1.0); + + materialStateEx.setDisplacement(0, 2.0); + materialStateEx.setDisplacement(1, 2.0); + materialStateEx.setDisplacement(2, 2.0); + materialStateEx.setRotation(0, 2.0); + materialStateEx.setRotation(1, 2.0); + materialStateEx.setRotation(2, 2.0); + materialStateEx.setSolidPressure(2.0); + + // Set pressure s.t. dP in TPSA source is nonzero + auto& sol = simulator->model().solution(/*timeIdx=*/0)[globI]; + sol[0] = 12.5e5; + simulator->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + + // Boundary info + // TODO: FIXED and FREE tests + BoundaryConditionData bcdata { Opm::BCMECHType::NONE, std::vector {0.0}, 0, 0.0 }; + + // + // Volume term + // + FVector volTerm; + LocalResidual::computeVolumeTerm(volTerm, + materialStateIn, + problem, + globI); + + for (std::size_t i = 0; i < volTerm.size(); ++i) { + if (i < 3) { + BOOST_CHECK_CLOSE(volTerm[i].value(), 0.0, 1.0e-3); + } + else if (i >= 3 && i < volTerm.size() - 1) { + // = rotation / SMODULUS + BOOST_CHECK_CLOSE(volTerm[i].value(), 1.0 / 3.0e9, 1.0e-3); + } + else { + // = solidPressure / LAME + BOOST_CHECK_CLOSE(volTerm[i].value(), 1.0 / 2.0e9, 1.0e-3); + } + } + + // + // Face term + // + FVector faceTerm; + LocalResidual::computeFaceTerm(faceTerm, + materialStateIn, + materialStateEx, + problem, + globI, + globJ); + + for (std::size_t i = 0; i < faceTerm.size(); ++i) { + if (i < 2) { + BOOST_CHECK_CLOSE(faceTerm[i].value(), -40000001.66666, 1.0e-3); + } + else if (i == 2) { + BOOST_CHECK_CLOSE(faceTerm[i].value(), -39999998.33333, 1.0e-3); + } + else if (i == 3) { + BOOST_CHECK_CLOSE(faceTerm[i].value(), 0.0, 1.0e-3); + } + else if (i == 4 || i == 6) { + BOOST_CHECK_CLOSE(faceTerm[i].value(), -1.33333, 1.0e-3); + } + else { + BOOST_CHECK_CLOSE(faceTerm[i].value(), 1.33333, 1.0e-3); + } + } + + // + // Source term + // + FVector srcTerm; + LocalResidual::computeSourceTerm(srcTerm, + problem, + globI, + /*timeIdx=*/0); + + for (std::size_t i = 0; i < srcTerm.size(); ++i) { + if (i == 6) { + BOOST_CHECK_CLOSE(srcTerm[i].value(), -0.000125, 1.0e-2); + } + else { + BOOST_CHECK_CLOSE(srcTerm[i].value(), 0.0, 1.0e-3); + } + } + + // + // Boundary term + // + FVector bndryTerm; + LocalResidual::computeBoundaryTerm(bndryTerm, + materialStateIn, + bcdata, + problem, + globI); + for (std::size_t i = 0; i < bndryTerm.size(); ++i) { + if (i < 2) { + BOOST_CHECK_CLOSE(bndryTerm[i].value(), 120000001.0, 1.0e-6); + } + else if (i == 2) { + BOOST_CHECK_CLOSE(bndryTerm[i].value(), 119999999.0, 1.0e-6); + } + else { + BOOST_CHECK_CLOSE(bndryTerm[i].value(), 0.0, 1.0e-6); + } + } +} diff --git a/tests/test_tpsa_primaryvariables.cpp b/tests/test_tpsa_primaryvariables.cpp new file mode 100644 index 0000000000..9df0383956 --- /dev/null +++ b/tests/test_tpsa_primaryvariables.cpp @@ -0,0 +1,81 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2025, NORCE AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include + +#define BOOST_TEST_MODULE TpsaPrimaryVariablesTests + +#include + +#include + +#include +#include + +#include +#include + +#include + + +namespace Opm::Properties::TTag { + struct TpsaTestTypeTag { + using InheritsFrom = std::tuple; + }; +} + +BOOST_AUTO_TEST_CASE(ElasticityPrimVarTest) { + using TypeTag = Opm::Properties::TTag::TpsaTestTypeTag; + using Scalar = Opm::GetPropType; + + // Check straightforward assignment + Opm::ElasticityPrimaryVariables priVars; + for (std::size_t i = 0; i < priVars.size(); ++i) { + priVars[i] = static_cast(i); + } + BOOST_CHECK_EQUAL(priVars.size(), 7U); + for (std::size_t i = 0; i < priVars.size(); ++i) { + const auto eval_var = priVars.makeEvaluation(i, 0); + BOOST_CHECK_CLOSE(eval_var.value(), static_cast(i), 1.0e-6); + BOOST_CHECK_CLOSE(eval_var.derivative(i), 1.0, 1.0e-6); + } + + // Assign from material state container + Opm::MaterialStateTPSA ms; + for (int i = 0; i < 3; ++i) { + ms.setDisplacement(i, static_cast(i) + 10.0); + ms.setRotation(i, static_cast(i) + 20.0); + } + ms.setSolidPressure(30.0); + Opm::ElasticityPrimaryVariables priVarsFromMs; + priVarsFromMs.assignNaive(ms); + for (std::size_t i = 0; i < 3; ++i) { + const auto eval_disp = priVarsFromMs.makeEvaluation(i, 0); + BOOST_CHECK_CLOSE(eval_disp.value(), ms.displacement(i), 1.0e-6); + BOOST_CHECK_CLOSE(eval_disp.derivative(i), 1.0, 1.0e-6); + + const auto eval_rot = priVarsFromMs.makeEvaluation(i + 3, 0); + BOOST_CHECK_CLOSE(eval_rot.value(), ms.rotation(i), 1.0e-6); + BOOST_CHECK_CLOSE(eval_rot.derivative(i + 3), 1.0, 1.0e-6); + } + const auto eval_sp = priVarsFromMs.makeEvaluation(6, 0); + BOOST_CHECK_CLOSE(eval_sp.value(), ms.solidPressure(), 1.0e-6); + BOOST_CHECK_CLOSE(eval_sp.derivative(6), 1.0, 1.0e-6); +} diff --git a/tests/tpsa_ex.data b/tests/tpsa_ex.data new file mode 100644 index 0000000000..e7f137804f --- /dev/null +++ b/tests/tpsa_ex.data @@ -0,0 +1,130 @@ +---------------- +RUNSPEC +---------------- + +TITLE + TWO CELLS X-DIR TPSA / + +DIMENS + 2 1 1 / + +EQLDIMS +/ + +TABDIMS +/ + +WELLDIMS +1 1 1 1 / + +WATER + +MECH + +MECHSOLV +TPSA / + +METRIC + +START + 1 JAN 2025 / + +UNIFIN +UNIFOUT + +---------------- +GRID +---------------- + +DX + 100 200 / +DY + 2*100 / +DZ + 2*100 / + +TOPS + 50 50 / + +PORO + 2*0.25 / +PERMX + 2*1000 / +COPY + PERMX PERMY / + PERMX PERMZ / +/ + +BIOTCOEF + 2*1.0 / + +LAME + 2*2.0 / + +SMODULUS + 2*3.0 / + +INIT + +---------------- +PROPS +---------------- + +ROCKOPTS + 1* STORE / + +ROCK + 1.0 1E-05 / + +DENSITY + 860.04 1033.0 0.853 / + +PVTW + 10 1.0 4E-5 0.5 / + +---------------- +SOLUTION +---------------- + +EQUIL + 100.0 10.0 / + +RTEMPVD + 100.0 25.0 + 200.0 25.0 / + +RPTRST + BASIC=2 / + +---------------- +SUMMARY +---------------- +FPR + +---------------- +SCHEDULE +---------------- + +RPTRST + BASIC=2 / + +WELSPECS + W1 INJ 1 1 1* WAT/ +/ + +COMPDAT + W1 1 1 1 1 OPEN 2* 0.2 / +/ + +WCONINJE + W1 WAT OPEN RATE 100.0 1* 900 / +/ + +--SOURCE +-- 1 1 1 WATER 100.0 / +--/ + +TSTEP + 2 / + +END