diff --git a/dumux/freeflow/navierstokes/model.hh b/dumux/freeflow/navierstokes/model.hh
index a4f6c02a9e64869d9fb1cf7bcbc4929917a17207..417604bddb29d945a44864ee461483b5e4fe9257 100644
--- a/dumux/freeflow/navierstokes/model.hh
+++ b/dumux/freeflow/navierstokes/model.hh
@@ -103,6 +103,12 @@ struct NavierStokesModelTraits
//! the indices
using Indices = NavierStokesIndices<dim()>;
+
+ //! return the name of the primary variables
+ static std::string primaryVariableName(int pvIdx, int state = 0)
+ {
+ return pvIdx == 0 ? "p" : "faceVelocity";
+ }
};
/*!
diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index b5f07746100b515253dde5a25cee1cbf5e903d09..aa2cdd60c03e0cd256780b877ccfadb3a6a52219 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -30,9 +30,11 @@
#include <unordered_set>
#include <dune/common/exceptions.hh>
+#include <dune/common/indices.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/typetraits/isvalid.hh>
+#include <dumux/common/typetraits/vector.hh>
#include <dumux/io/vtk/vtkreader.hh>
namespace Dumux {
@@ -119,6 +121,65 @@ auto loadSolutionFromVtkFile(const std::string fileName,
}
}
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to read from a file into a solution vector
+ */
+template <class SolutionVector, class PvNamesFunc>
+auto loadSolutionFromVtkFile(const std::string fileName,
+ const VTKReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
+-> typename std::enable_if_t<decltype(isMultiTypeBlockVector<SolutionVector>())::value, void>
+{}
+
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to read from two files into a staggered solution vector
+ */
+template <class SolutionVector, class PvNamesFunc>
+auto loadStaggeredSolutionFromVtkFiles(const std::string fileName,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
+-> typename std::enable_if_t<decltype(!isMultiTypeBlockVector<SolutionVector>())::value, void>
+{}
+
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to read from two files into a staggered solution vector
+ */
+template <class SolutionVector, class PvNamesFunc>
+auto loadStaggeredSolutionFromVtkFiles(const std::string baseFileName,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
+-> typename std::enable_if_t<decltype(isMultiTypeBlockVector<SolutionVector>())::value, void>
+{
+
+ // assume that the first component contains the cell data
+ auto& cellSol = sol[Dune::index_constant<0>{}];
+ using CellPrimaryVariables = typename std::decay_t<decltype(cellSol)>::block_type;
+ using Scalar = typename CellPrimaryVariables::field_type;
+ VTKReader cellVtk(baseFileName + ".vtu");
+ for (size_t pvIdx = 0; pvIdx < CellPrimaryVariables::dimension; ++pvIdx)
+ {
+ const auto vec = cellVtk.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx),
+ VTKReader::DataType::cellData);
+ for (size_t i = 0; i < cellSol.size(); ++i)
+ cellSol[i][pvIdx] = vec[i];
+ }
+
+ // assume that the second component contains the face data
+ auto& faceSol = sol[Dune::index_constant<1>{}];
+ using FacePrimaryVariables = typename std::decay_t<decltype(faceSol)>::block_type;
+ auto nameSize = baseFileName.size();
+ // assume that baseFileName contains numbers like '-000123' in the end
+ VTKReader faceVtk(baseFileName.substr(0, nameSize - 6) + "-face" + baseFileName.substr(nameSize - 6) + ".vtp");
+ const auto vec = faceVtk.readData<std::vector<Scalar>>(pvNamesFunc(CellPrimaryVariables::dimension),
+ VTKReader::DataType::pointData);
+ for (size_t i = 0; i < faceSol.size(); ++i)
+ faceSol[i][0] = std::accumulate(&vec[3*i], &vec[3*(i+1)], 0.0);
+}
+
/*!
* \ingroup InputOutput
* \brief helper function to determine the primray variable names of a model with privar state
@@ -129,12 +190,12 @@ std::string primaryVariableName(int pvIdx, int state)
static auto numStates = (1 << ModelTraits::numPhases()) - 1;
const auto paramNameWithState = "LoadSolution.PriVarNamesState" + std::to_string(state);
- if (haveParam("LoadSolution.PriVarNames") && !haveParam(paramNameWithState))
+ if (hasParam("LoadSolution.PriVarNames") && !hasParam(paramNameWithState))
{
DUNE_THROW(Dune::NotImplemented, "please provide LoadSolution.PriVarNamesState1..." << numStates
<< " or remove LoadSolution.PriVarNames to use default names");
}
- else if (haveParam(paramNameWithState))
+ else if (hasParam(paramNameWithState))
{
const auto pvNames = getParam<std::vector<std::string>>(paramNameWithState);
return pvNames[pvIdx];
@@ -150,7 +211,7 @@ std::string primaryVariableName(int pvIdx, int state)
template<class ModelTraits>
std::string primaryVariableName(int pvIdx)
{
- if (haveParam("LoadSolution.PriVarNames"))
+ if (hasParam("LoadSolution.PriVarNames"))
{
static auto pvNames = getParam<std::vector<std::string>>("LoadSolution.PriVarNames");
return pvNames[pvIdx];
@@ -179,6 +240,8 @@ void loadSolution(const std::string& fileName,
? VTKReader::DataType::pointData : VTKReader::DataType::cellData;
loadSolutionFromVtkFile(fileName, dataType, pvNamesFunc, sol);
}
+ else if (extension == fileName && discMethod == DiscretizationMethod::staggered)
+ loadStaggeredSolutionFromVtkFiles(fileName, pvNamesFunc, sol);
else
DUNE_THROW(Dune::NotImplemented, "loadSolution for extension " << extension);
}
diff --git a/test/freeflow/navierstokes/CMakeLists.txt b/test/freeflow/navierstokes/CMakeLists.txt
index d74193f31913083ef9ac6272c0f6b4870c41e001..cdcf5b0a93347c16868bda28471c476ae417aad9 100644
--- a/test/freeflow/navierstokes/CMakeLists.txt
+++ b/test/freeflow/navierstokes/CMakeLists.txt
@@ -70,7 +70,20 @@ dune_add_test(NAME test_channel_navierstokes
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/channel-navierstokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes-00002.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes")
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes -Vtk.WriteFaceData 1")
+
+dune_add_test(NAME test_channel_navierstokes_restart
+ SOURCES test_channel.cc
+ COMPILE_DEFINITIONS ENABLE_NAVIERSTOKES=1
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/channel-navierstokes-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes_restart-00001.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes -Vtk.WriteFaceData 1 -TimeLoop.DtInitial 1 -Restart.Time 1 -Restart.File test_channel_navierstokes-00001 -Problem.Name test_channel_navierstokes_restart")
+
+# the restart test has to run after the test that produces the corresponding vtu file
+set_tests_properties(test_channel_navierstokes_restart PROPERTIES DEPENDS test_channel_navierstokes)
dune_add_test(NAME test_stokes_donea_no_caching
SOURCES test_donea.cc
diff --git a/test/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc
index bee7ed183fd72be3295274245c615602bc303e85..57d5345763ac9c51745d83b45aebe62ea1951220 100644
--- a/test/freeflow/navierstokes/test_channel.cc
+++ b/test/freeflow/navierstokes/test_channel.cc
@@ -50,6 +50,7 @@
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
+#include <dumux/io/loadsolution.hh>
#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
@@ -129,14 +130,7 @@ int main(int argc, char** argv) try
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");
-
- // instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
- problem->setTimeLoop(timeLoop);
+ Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -145,9 +139,21 @@ int main(int argc, char** argv) try
SolutionVector x;
x[FVGridGeometry::cellCenterIdx()].resize(numDofsCellCenter);
x[FVGridGeometry::faceIdx()].resize(numDofsFace);
- problem->applyInitialSolution(x);
+ if (restartTime > 0)
+ {
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ auto fileName = getParam<std::string>("Restart.File");
+ loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits>, x);
+ }
+ else
+ problem->applyInitialSolution(x);
auto xOld = x;
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ problem->setTimeLoop(timeLoop);
+
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
@@ -157,7 +163,7 @@ int main(int argc, char** argv) try
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
- vtkWriter.write(0.0);
+ vtkWriter.write(restartTime);
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;