TPSA implementation with coupling to Flow.

- TPSA-Flow coupling done at the nonlinear solver step level:
- -- Lagged coupling -> TPSA lags behind one time step
- -- Fixed-stress -> Iterate until Flow and TPSA converge
- Generic TPSA solver using Newton with associated linearizer.
- Linear elasticity equations implemented
- Specialized version of ISTL linear solver for TPSA
- Output to restart files and VTK
This commit is contained in:
Svenn Tveit committed 2026-03-26 15:23:47 +01:00
1 parent 2f7d57b7b3
commit 106d8d51c9
51 files changed
+6964 -33

No files matched your search

+192
View File
@@ -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 <http://www.gnu.org/licenses/>.
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>
#define BOOST_TEST_MODULE TestFacePropertiesTPSA
#define BOOST_TEST_NO_MAIN
#include <boost/test/unit_test.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/input/eclipse/Parser/Parser.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/flow/FacePropertiesTPSA.hpp>
#include <string>
#include <vector>
#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<char> 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<GridView>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
using FacePropertiesTPSA = Opm::FacePropertiesTPSA<Grid, GridView, ElementMapper, CartesianIndexMapper, double>;
using DimVector = Dune::FieldVector<double, GridView::dimensionworld>;
// 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>(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<DimVector> 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);
}
+284
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#define BOOST_TEST_MODULE TpsaLocalResidualTests
#include <boost/test/unit_test.hpp>
#include <dune/common/fvector.hh>
#include <opm/input/eclipse/Schedule/BCProp.hpp>
#include <opm/material/materialstates/MaterialStateTPSA.hpp>
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/BlackoilModelProperties.hpp>
#include <opm/simulators/flow/FlowProblemTPSA.hpp>
#include <opm/simulators/flow/TTagFlowProblemTPSA.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
#include <tuple>
#include <vector>
namespace Opm::Properties {
namespace TTag {
struct TpsaTestTypeTag {
using InheritsFrom = std::tuple<FlowProblem, FlowProblemTpsa>;
};
}
template <class TypeTag>
struct Problem<TypeTag, TTag::TpsaTestTypeTag>
{ using type = FlowProblemTPSA<TypeTag>; };
template <class TypeTag>
struct EnableMech<TypeTag, TTag::TpsaTestTypeTag>
{ static constexpr bool value = true; };
template<class TypeTag>
struct Indices<TypeTag, TTag::TpsaTestTypeTag>
{
private:
using BaseTypeTag = TTag::FlowProblem;
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
public:
using type = BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
/*enabledCompIdx=*/FluidSystem::waterCompIdx,
getPropValue<TypeTag, Properties::EnableBioeffects>()>;
};
}
namespace {
template <class TypeTag>
std::unique_ptr<Opm::GetPropType<TypeTag, Opm::Properties::Simulator>>
initSimulator(const char *filename)
{
using Simulator = Opm::GetPropType<TypeTag, Opm::Properties::Simulator>;
const auto filenameArg = std::string {"--ecl-deck-file-name="} + filename;
const char* argv[] = {
"test_tpsa",
filenameArg.c_str()
};
Opm::Parameters::reset();
Opm::registerAllParameters_<TypeTag>(false);
Opm::BlackoilModelParameters<double>::registerParameters();
Opm::registerEclTimeSteppingParameters<double>();
Opm::Parameters::Register<Opm::Parameters::EnableTerminalOutput>("Do *NOT* use!");
Opm::Parameters::SetDefault<Opm::Parameters::ThreadsPerProcess>(1);
Opm::Parameters::endRegistration();
Opm::setupParameters_<TypeTag>(/*argc=*/sizeof(argv) / sizeof(argv[0]),
argv,
/*registerParams=*/false,
/*allowUnused*/false,
/*handleHelp*/true,
/*myRank*/0);
Opm::FlowGenericVanguard::readDeck(filename);
return std::make_unique<Simulator>();
}
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<Opm::Parallel::Communication>());
}
};
struct BoundaryConditionData
{
Opm::BCMECHType type;
std::vector<double> 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<TypeTag>;
using Evaluation = Opm::GetPropType<TypeTag, Opm::Properties::EvaluationTPSA>;
using MaterialState = Opm::MaterialStateTPSA<Evaluation>;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
using FVector = Dune::FieldVector<Evaluation, 7>;
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
// Setup
auto simulator = initSimulator<TypeTag>("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<double> {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);
}
}
}
+81
View File
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#define BOOST_TEST_MODULE TpsaPrimaryVariablesTests
#include <boost/test/unit_test.hpp>
#include <opm/material/materialstates/MaterialStateTPSA.hpp>
#include <opm/models/tpsa/elasticityprimaryvariables.hpp>
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/flow/BlackoilModelProperties.hpp>
#include <opm/simulators/flow/TTagFlowProblemTPSA.hpp>
#include <tuple>
namespace Opm::Properties::TTag {
struct TpsaTestTypeTag {
using InheritsFrom = std::tuple<FlowProblem, FlowProblemTpsa>;
};
}
BOOST_AUTO_TEST_CASE(ElasticityPrimVarTest) {
using TypeTag = Opm::Properties::TTag::TpsaTestTypeTag;
using Scalar = Opm::GetPropType<TypeTag, Opm::Properties::Scalar>;
// Check straightforward assignment
Opm::ElasticityPrimaryVariables<TypeTag> priVars;
for (std::size_t i = 0; i < priVars.size(); ++i) {
priVars[i] = static_cast<double>(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<double>(i), 1.0e-6);
BOOST_CHECK_CLOSE(eval_var.derivative(i), 1.0, 1.0e-6);
}
// Assign from material state container
Opm::MaterialStateTPSA<Scalar> ms;
for (int i = 0; i < 3; ++i) {
ms.setDisplacement(i, static_cast<double>(i) + 10.0);
ms.setRotation(i, static_cast<double>(i) + 20.0);
}
ms.setSolidPressure(30.0);
Opm::ElasticityPrimaryVariables<TypeTag> 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);
}
+130
View File
@@ -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