From d82bca6e6dc37b15e5e3997c5b7b1697594c2e44 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 26 Jul 2018 13:07:01 +0200
Subject: [PATCH] [io][vtk][parallel] read from parallel Vtk files
Works for non-staggered discretizations and solution vectors without
state so far.
---
dumux/io/loadsolution.hh | 182 ++++++++++++++++--
dumux/io/vtk/vtkreader.hh | 56 +++++-
test/freeflow/navierstokes/test_channel.cc | 4 +-
.../2p/implicit/incompressible/test_2p_fv.cc | 2 +-
.../2p2c/implicit/test_2p2c_fv.cc | 2 +-
5 files changed, 224 insertions(+), 22 deletions(-)
diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index dad597320e..c047c386d3 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -51,15 +51,61 @@ struct hasState
} // end namespace Detail
+template <class FieldType, class Container, class GridView, int codim>
+class LoadSolutionDataHandle
+: public Dune::CommDataHandleIF< LoadSolutionDataHandle<FieldType, Container, GridView, codim>,
+ FieldType >
+{
+ using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+
+public:
+ LoadSolutionDataHandle(Container &container,
+ const GridView &gridView)
+ : mapper_(gridView, Dune::mcmgLayout(Dune::Codim<codim>()))
+ , container_(container)
+ {}
+
+ bool contains(int dim, int cd) const
+ { return cd == codim; }
+
+ bool fixedsize(int dim, int cd) const
+ { return true; }
+
+ template<class EntityType>
+ size_t size (const EntityType &e) const
+ { return 1; }
+
+ template<class MessageBufferImp, class EntityType>
+ void gather(MessageBufferImp &buff, const EntityType &e) const
+ {
+ int vIdx = mapper_.index(e);
+ buff.write(container_[vIdx]);
+ }
+
+ template<class MessageBufferImp, class EntityType>
+ void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
+ {
+ int vIdx = mapper_.index(e);
+
+ FieldType tmp;
+ buff.read(tmp);
+ container_[vIdx] = tmp;
+ }
+
+private:
+ ElementMapper mapper_;
+ Container &container_;
+};
+
/*!
* \ingroup InputOutput
- * \brief helper function to read from a file into a solution vector
+ * \brief read from a sequential file into a solution vector without state
*/
template <class SolutionVector, class PvNamesFunc>
-auto loadSolutionFromVtkFile(const std::string fileName,
- const VTKReader::DataType& dataType,
- PvNamesFunc&& pvNamesFunc,
- SolutionVector& sol)
+auto loadSolutionFromSequentialVtkFile(const std::string fileName,
+ const VTKReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
-> typename std::enable_if_t<!decltype(isValid(Detail::hasState())(sol[0]))::value, void>
{
VTKReader vtu(fileName);
@@ -80,13 +126,73 @@ auto loadSolutionFromVtkFile(const std::string fileName,
/*!
* \ingroup InputOutput
- * \brief helper function to read from a file into a solution vector
+ * \brief read from a parallel file into a solution vector without state
+ */
+template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
+auto loadSolutionFromParallelVtkFile(const std::string fileName,
+ const VTKReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol,
+ const FVGridGeometry& fvGridGeometry)
+-> typename std::enable_if_t<!decltype(isValid(Detail::hasState())(sol[0]))::value, void>
+{
+ VTKReader vtu(fileName);
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using Scalar = typename PrimaryVariables::field_type;
+ using GridView = typename FVGridGeometry::GridView;
+ static constexpr auto dim = GridView::dimension;
+
+ for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+ {
+ const auto pvName = pvNamesFunc(pvIdx);
+ auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
+ if (vec.size() != sol.size() && fvGridGeometry.gridView().comm().size() == 1)
+ {
+ DUNE_THROW(Dune::IOError, "Size mismatch between solution vector and read data ("
+ << sol.size() << " != " << vec.size() << ")");
+ }
+
+ std::vector<bool> visited;
+ if (dataType == VTKReader::DataType::pointData)
+ {
+ visited.resize(fvGridGeometry.gridView().size(dim), false);
+ }
+
+ std::size_t i = 0;
+ for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
+ {
+ if (dataType == VTKReader::DataType::cellData)
+ {
+ auto eIdx = fvGridGeometry.elementMapper().index(element);
+ sol[eIdx][pvIdx] = vec[i];
+ ++i;
+ }
+ else
+ {
+ for (int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
+ {
+ auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
+ if (!visited[vIdxGlobal])
+ {
+ sol[vIdxGlobal][pvIdx] = vec[i];
+ ++i;
+ visited[vIdxGlobal] = true;
+ }
+ }
+ }
+ }
+ }
+}
+
+/*!
+ * \ingroup InputOutput
+ * \brief read from a sequential file into a solution vector with state
*/
template <class SolutionVector, class PvNamesFunc>
-auto loadSolutionFromVtkFile(const std::string fileName,
- const VTKReader::DataType& dataType,
- PvNamesFunc&& pvNamesFunc,
- SolutionVector& sol)
+auto loadSolutionFromSequentialVtkFile(const std::string fileName,
+ const VTKReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
-> typename std::enable_if_t<decltype(isValid(Detail::hasState())(sol[0]))::value, void>
{
VTKReader vtu(fileName);
@@ -121,6 +227,21 @@ auto loadSolutionFromVtkFile(const std::string fileName,
}
}
+/*!
+ * \ingroup InputOutput
+ * \brief read from a parallel file into a solution vector with state
+ */
+template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
+auto loadSolutionFromParallelVtkFile(const std::string fileName,
+ const VTKReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol,
+ const FVGridGeometry& fvGridGeometry)
+-> typename std::enable_if_t<decltype(isValid(Detail::hasState())(sol[0]))::value, void>
+{
+ DUNE_THROW(Dune::NotImplemented, "reading solution with state from a parallel Vtk file");
+}
+
/*!
* \ingroup InputOutput
* \brief helper function to determine the primary variable names of a model with privar state
@@ -223,25 +344,56 @@ std::string primaryVariableNameFace(int pvIdx)
* \brief load a solution vector from file
* \note Supports the following file extensions: *.vtu *.vtp
*/
-template <class SolutionVector, class PvNamesFunc>
+template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
void loadSolution(const std::string& fileName,
DiscretizationMethod discMethod,
PvNamesFunc&& pvNamesFunc,
- SolutionVector& sol)
+ SolutionVector& sol,
+ const FVGridGeometry& fvGridGeometry)
{
const auto extension = fileName.substr(fileName.find_last_of(".") + 1);
+ auto dataType = discMethod == DiscretizationMethod::box ?
+ VTKReader::DataType::pointData : VTKReader::DataType::cellData;
if (extension == "vtu" || extension == "vtp")
{
- auto dataType = discMethod == DiscretizationMethod::box ?
- VTKReader::DataType::pointData : VTKReader::DataType::cellData;
if (discMethod == DiscretizationMethod::staggered && extension == "vtp")
dataType = VTKReader::DataType::pointData;
- loadSolutionFromVtkFile(fileName, dataType, pvNamesFunc, sol);
+ loadSolutionFromSequentialVtkFile(fileName, dataType, pvNamesFunc, sol);
+ }
+ else if (extension == "pvtu" || extension == "pvtp")
+ {
+ if (discMethod == DiscretizationMethod::staggered)
+ DUNE_THROW(Dune::NotImplemented, "reading staggered solution from a parallel Vtk file");
+
+ loadSolutionFromParallelVtkFile(fileName, dataType, pvNamesFunc, sol, fvGridGeometry);
}
else
DUNE_THROW(Dune::NotImplemented, "loadSolution for extension " << extension);
+
+ // synchronize values on ghost and overlap dofs
+ if (fvGridGeometry.gridView().comm().size() > 1)
+ {
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using GridView = typename FVGridGeometry::GridView;
+ if (dataType == VTKReader::DataType::cellData)
+ {
+ LoadSolutionDataHandle<PrimaryVariables, SolutionVector, GridView, 0>
+ dataHandle(sol, fvGridGeometry.gridView());
+ fvGridGeometry.gridView().communicate(dataHandle,
+ Dune::InteriorBorder_All_Interface,
+ Dune::ForwardCommunication);
+ }
+ else
+ {
+ LoadSolutionDataHandle<PrimaryVariables, SolutionVector, GridView, GridView::dimension>
+ dataHandle(sol, fvGridGeometry.gridView());
+ fvGridGeometry.gridView().communicate(dataHandle,
+ Dune::InteriorBorder_All_Interface,
+ Dune::ForwardCommunication);
+ }
+ }
}
} // namespace Dumux
diff --git a/dumux/io/vtk/vtkreader.hh b/dumux/io/vtk/vtkreader.hh
index 30ccb3fa28..127ce98b80 100644
--- a/dumux/io/vtk/vtkreader.hh
+++ b/dumux/io/vtk/vtkreader.hh
@@ -31,6 +31,7 @@
#include <type_traits>
#include <unordered_map>
+#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/io/file/vtk/common.hh>
@@ -60,10 +61,13 @@ public:
explicit VTKReader(const std::string& fileName)
: fileName_(fileName)
{
+ if (Dune::MPIHelper::getCollectiveCommunication().size() > 1)
+ getSequentialNameFromParallelName_();
+
using namespace tinyxml2;
- const auto eResult = doc_.LoadFile(fileName.c_str());
+ const auto eResult = doc_.LoadFile(fileName_.c_str());
if (eResult != XML_SUCCESS)
- DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
+ DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName_ << ".");
}
/*!
@@ -133,6 +137,52 @@ public:
}
private:
+ /*!
+ * \brief get sequential Vtk filename from a parallel one
+ */
+ void getSequentialNameFromParallelName_()
+ {
+ using namespace tinyxml2;
+
+ const auto eResult = doc_.LoadFile(fileName_.c_str());
+ if (eResult != XML_SUCCESS)
+ DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName_ << ".");
+
+ const XMLElement* pieceNode = doc_.FirstChildElement("VTKFile");
+ if (pieceNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get 'VTKFile' node in " << fileName_ << ".");
+
+ pieceNode = pieceNode->FirstChildElement("PUnstructuredGrid");
+ if (pieceNode == nullptr)
+ pieceNode = doc_.FirstChildElement("VTKFile")->FirstChildElement("PPolyData");
+ if (pieceNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get 'PUnstructuredGrid' or 'PPolyData' node in " << fileName_ << ".");
+
+ pieceNode = pieceNode->FirstChildElement("Piece");
+ const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+ for (int rank = 0; rank < comm.size(); ++rank)
+ {
+ if (pieceNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get 'Piece' node for rank "
+ << rank << " in " << fileName_ << ".");
+
+ const char *seqName = pieceNode->Attribute("Source");
+
+ if (seqName == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get 'Source' attribute of a 'Piece' node in " << fileName_);
+
+ if (comm.rank() == rank)
+ {
+ fileName_ = seqName;
+ break;
+ }
+
+ pieceNode = pieceNode->NextSiblingElement("Piece");
+ }
+
+ doc_.Clear();
+ }
+
/*!
* \brief Read a grid from a vtk/vtu/vtp file
* \param factory the (emtpy) grid factory
@@ -387,7 +437,7 @@ private:
adaptPointDimension_(std::vector<Dune::FieldVector<double, 3>>&& points3D) const
{ return points3D; }
- const std::string fileName_; //!< the vtk file name
+ std::string fileName_; //!< the vtk file name
tinyxml2::XMLDocument doc_; //!< the xml document created from file with name fileName_
};
diff --git a/test/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc
index 7350a12d50..9b85acf622 100644
--- a/test/freeflow/navierstokes/test_channel.cc
+++ b/test/freeflow/navierstokes/test_channel.cc
@@ -145,11 +145,11 @@ int main(int argc, char** argv) try
auto fileNameCell = getParam<std::string>("Restart.FileCell");
loadSolution(fileNameCell, FVGridGeometry::discMethod,
- primaryVariableNameCell<ModelTraits>, x[Dune::index_constant<0>{}]);
+ primaryVariableNameCell<ModelTraits>, x[Dune::index_constant<0>{}], *fvGridGeometry);
auto fileNameFace = getParam<std::string>("Restart.FileFace");
loadSolution(fileNameFace, FVGridGeometry::discMethod,
- primaryVariableNameFace<ModelTraits>, x[Dune::index_constant<1>{}]);
+ primaryVariableNameFace<ModelTraits>, x[Dune::index_constant<1>{}], *fvGridGeometry);
}
else
problem->applyInitialSolution(x);
diff --git a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
index 6098c4a07e..06f13ebd40 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
+++ b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
@@ -138,7 +138,7 @@ int main(int argc, char** argv) try
{
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
auto fileName = getParam<std::string>("Restart.File");
- loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits>, x);
+ loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits>, x, *fvGridGeometry);
}
else
problem->applyInitialSolution(x);
diff --git a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
index 3ef717512c..a1190fff94 100644
--- a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
+++ b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
@@ -107,7 +107,7 @@ int main(int argc, char** argv) try
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
auto fileName = getParam<std::string>("Restart.File");
- loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits, FluidSystem>, x);
+ loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits, FluidSystem>, x, *fvGridGeometry);
}
else
problem->applyInitialSolution(x);
--
GitLab