From 79601d68b2537ffb2a5da4bf37bf9d0118a061ca Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 4 Jul 2018 10:31:51 +0200
Subject: [PATCH] [io][restart] make restart work for 2p2c
Distinguish between switchable primary variables and standard ones.
Outsource filling of the solution vector to a correspondingly
specialized function `setPrimaryVariables`.
---
dumux/io/restart.hh | 136 ++++++++++++++----
.../2p2c/implicit/CMakeLists.txt | 9 ++
.../2p2c/implicit/injectionproblem.hh | 8 ++
.../2p2c/implicit/test_2p2c_fv.cc | 24 ++--
4 files changed, 140 insertions(+), 37 deletions(-)
diff --git a/dumux/io/restart.hh b/dumux/io/restart.hh
index 5a1103d33a..9848135863 100644
--- a/dumux/io/restart.hh
+++ b/dumux/io/restart.hh
@@ -35,9 +35,20 @@
#include <fstream>
#include <sstream>
+#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/io/tinyxml2/tinyxml2.h>
namespace Dumux {
+
+//! helper struct detecting if a PrimaryVariables object has a state() function
+struct hasState
+{
+ template<class PrimaryVariables>
+ auto operator()(PrimaryVariables&& priVars)
+ -> decltype(priVars.state())
+ {}
+};
+
/*!
* \ingroup InputOutput
* \brief Load or save a state of a model to/from the harddisk.
@@ -289,10 +300,36 @@ public:
"Restart::restartFileList()");
}
+ template <class FVGridGeometry, class SolutionVector>
+ static void loadSolutionFromVtkFile(const FVGridGeometry& fvGridGeometry,
+ const std::vector<std::string>& pvNames,
+ SolutionVector& sol)
+ {
+ using namespace tinyxml2;
+ auto fileName = getParam<std::string>("Restart.File");
+
+ XMLDocument xmlDoc;
+ auto eResult = xmlDoc.LoadFile(fileName.c_str());
+ if (eResult != XML_SUCCESS)
+ DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
+
+ XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
+ if (pieceNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName << ".");
+
+ XMLElement *cellDataNode = pieceNode->FirstChildElement("CellData");
+ if (cellDataNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get CellData node in " << fileName << ".");
+
+ setPrimaryVariables_(cellDataNode, fvGridGeometry, pvNames, sol);
+ }
+
+private:
+
template <class Scalar, class FVGridGeometry>
- static std::vector<Scalar> extractDataToVector(tinyxml2::XMLElement *xmlNode,
- const std::string& name,
- const FVGridGeometry& fvGridGeometry)
+ static std::vector<Scalar> extractDataToVector_(tinyxml2::XMLElement *xmlNode,
+ const std::string& name,
+ const FVGridGeometry& fvGridGeometry)
{
// loop over XML node siblings to find the correct data array
tinyxml2::XMLElement *dataArray = xmlNode->FirstChildElement("DataArray");
@@ -319,39 +356,90 @@ public:
return vec;
}
- template <class FVGridGeometry, class SolutionVector>
- static void loadSolutionFromVtkFile(const FVGridGeometry& fvGridGeometry,
- std::vector<std::string> pvNames,
- SolutionVector& sol)
+ template<class SolutionVector, class FVGridGeometry>
+ static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
+ const FVGridGeometry& fvGridGeometry,
+ std::vector<std::string> pvNames,
+ SolutionVector& sol)
+ -> typename std::enable_if_t<!decltype(isValid(hasState())(sol[0]))::value, void>
{
- using namespace tinyxml2;
- auto fileName = getParam<std::string>("Restart.File");
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using Scalar = typename PrimaryVariables::field_type;
- XMLDocument xmlDoc;
- auto eResult = xmlDoc.LoadFile(fileName.c_str());
- if (eResult != XML_SUCCESS)
- DUNE_THROW(Dune::IOError, "Couldn't open XML file.");
+ for (size_t pvIdx = 0; pvIdx < pvNames.size(); ++pvIdx)
+ {
+ auto pvName = pvNames[pvIdx];
- XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
- if (pieceNode == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't get Piece node.");
+ auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
- XMLElement *cellDataNode = pieceNode->FirstChildElement("CellData");
- if (cellDataNode == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't get CellData node.");
+ for (size_t i = 0; i < sol.size(); ++i)
+ sol[i][pvIdx] = vec[i];
+ }
+ }
+
+ template<class SolutionVector, class FVGridGeometry>
+ static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
+ const FVGridGeometry& fvGridGeometry,
+ std::vector<std::string> pvNames,
+ SolutionVector& sol)
+ -> typename std::enable_if_t<decltype(isValid(hasState())(sol[0]))::value, void>
+ {
+ auto phasePresenceIt = std::find(pvNames.begin(), pvNames.end(), "phase presence");
+ if (phasePresenceIt != pvNames.end())
+ {
+ auto vec = extractDataToVector_<int>(node, "phase presence", fvGridGeometry);
+
+ for (size_t i = 0; i < sol.size(); ++i)
+ sol[i].setState(vec[i]);
+
+ pvNames.erase(phasePresenceIt);
+ }
using PrimaryVariables = typename SolutionVector::block_type;
using Scalar = typename PrimaryVariables::field_type;
- for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+ for (size_t pvIdx = 0; pvIdx < pvNames.size(); ++pvIdx)
{
- auto vec = extractDataToVector<Scalar>(cellDataNode, pvNames[pvIdx], fvGridGeometry);
+ auto pvName = pvNames[pvIdx];
- for (size_t i = 0; i < sol.size(); ++i)
- sol[i][pvIdx] = vec[i];
+ auto switchedPvNames = getSwitchedPvNames_(pvName);
+
+ if (switchedPvNames[0] == pvName)
+ {
+ auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
+
+ for (size_t i = 0; i < sol.size(); ++i)
+ sol[i][pvIdx] = vec[i];
+ }
+ else
+ {
+ std::vector<std::vector<Scalar>> switchedPvsSol;
+ for (size_t switchedPvIdx = 0; switchedPvIdx < switchedPvNames.size(); ++switchedPvIdx)
+ switchedPvsSol.push_back(extractDataToVector_<Scalar>(node,
+ switchedPvNames[switchedPvIdx],
+ fvGridGeometry));
+
+ for (size_t i = 0; i < sol.size(); ++i)
+ sol[i][pvIdx] = switchedPvsSol[sol[i].state()-1][i];
+ }
}
}
-private:
+ static std::vector<std::string> getSwitchedPvNames_(std::string pvName)
+ {
+ std::vector<std::string> switchedPvNames;
+ int lastPos = -1;
+ size_t pos = 0;
+ do
+ {
+ pos = pvName.find('/', lastPos + 1);
+ switchedPvNames.push_back(pvName.substr(lastPos+1, pos - lastPos - 1));
+ std::cout << "Found switched PV name " << pvName.substr(lastPos+1, pos - lastPos - 1) << std::endl;
+ lastPos = pos;
+ } while (pos != std::string::npos);
+
+ return switchedPvNames;
+ }
+
std::string fileName_;
std::ifstream inStream_;
std::ofstream outStream_;
diff --git a/test/porousmediumflow/2p2c/implicit/CMakeLists.txt b/test/porousmediumflow/2p2c/implicit/CMakeLists.txt
index 2ccc8da3c9..7e5a51728e 100644
--- a/test/porousmediumflow/2p2c/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2p2c/implicit/CMakeLists.txt
@@ -20,6 +20,15 @@ dune_add_test(NAME test_2p2c_tpfa
${CMAKE_CURRENT_BINARY_DIR}/injection_tpfa-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p2c_tpfa test_2p2c_fv.input -Problem.Name injection_tpfa")
+dune_add_test(NAME test_2p2c_tpfa_restart
+ SOURCES test_2p2c_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=InjectionCCTpfaTypeTag ENABLECACHING=0
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/injectioncc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/injection_tpfa_restart-00004.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p2c_tpfa test_2p2c_fv.input -Problem.Name injection_tpfa_restart -TimeLoop.DtInitial 1488.5 -Restart.Time 2158.85 -Restart.File injection_tpfa-00004.vtu")
+
dune_add_test(NAME test_2p2c_mpfa
SOURCES test_2p2c_fv.cc
COMPILE_DEFINITIONS TYPETAG=InjectionCCMpfaTypeTag ENABLECACHING=0
diff --git a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
index 56e2d8093c..b23b3e03c8 100644
--- a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
@@ -30,6 +30,8 @@
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/box/properties.hh>
+#include <dumux/io/restart.hh>
+
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/2p2c/model.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
@@ -292,6 +294,12 @@ public:
// \}
+ template <class SolutionVector>
+ void applyRestartSolution(SolutionVector& sol) const
+ {
+ Restart::loadSolutionFromVtkFile(this->fvGridGeometry(), {"pw", "x_w^N2/x_n^H2O/Sn", "phase presence"}, sol);
+ }
+
private:
/*!
* \brief Evaluates the initial values for a control volume
diff --git a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
index efef6227db..1e80cd70e3 100644
--- a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
+++ b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
@@ -89,10 +89,19 @@ int main(int argc, char** argv) try
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
+
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(fvGridGeometry->numDofs());
- problem->applyInitialSolution(x);
+ problem->applyInitialSolution(x, restartTime);
auto xOld = x;
// the grid variables
@@ -100,24 +109,13 @@ int main(int argc, char** argv) try
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
- // get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- auto dt = getParam<Scalar>("TimeLoop.DtInitial");
-
- // check if we are about to restart a previously interrupted simulation
- Scalar restartTime = 0;
- if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
- restartTime = getParam<Scalar>("TimeLoop.Restart");
-
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
- vtkWriter.write(0.0);
+ vtkWriter.write(restartTime);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
--
GitLab