From 0e4b2e3c7bd723ba9be9016423801a6da1e33dae Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Mon, 19 Dec 2016 15:42:36 +0100
Subject: [PATCH] [vtk] Introduce new output module featuring standardized
output fields
---
dumux/implicit/problem.hh | 20 ++-
dumux/io/vtkoutputmodule.hh | 349 ++++++++++++++++++++++++++++++++++++
2 files changed, 368 insertions(+), 1 deletion(-)
create mode 100644 dumux/io/vtkoutputmodule.hh
diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh
index 24675fb4af..b10783c782 100644
--- a/dumux/implicit/problem.hh
+++ b/dumux/implicit/problem.hh
@@ -27,6 +27,8 @@
#include "model.hh"
#include <dumux/io/restart.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+
#include <dumux/implicit/adaptive/gridadapt.hh>
#include <dumux/common/boundingboxtree.hh>
@@ -138,6 +140,8 @@ public:
*/
void init()
{
+ vtkOutputModule_ = std::make_shared<VtkOutputModule<TypeTag>>(asImp_());
+
// set the initial condition of the model
model().init(asImp_());
@@ -901,6 +905,12 @@ public:
void addOutputVtkFields()
{}
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+ void addVtkOutputFields(VtkOutputModule<TypeTag>& outputModule) const
+ {}
+
/*!
* \brief Write the relevant secondary variables of the current
* solution into an VTK output file.
@@ -918,6 +928,8 @@ public:
model().addOutputVtkFields(model().curSol(), *resultWriter_);
asImp_().addOutputVtkFields();
resultWriter_->endWrite();
+
+ vtkOutputModule_->write(t);
}
/*!
@@ -1035,6 +1047,11 @@ public:
void postAdapt()
{}
+ VtkOutputModule<TypeTag>& vtkOutputModule() const
+ {
+ return *vtkOutputModule_;
+ }
+
protected:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
@@ -1057,7 +1074,7 @@ protected:
return *resultWriter_;
}
- //! Compute the point source map, i.e. which scvs have point source contributions
+ //! Compute the point source map, i.e. which scvs have point source contributions
void computePointSourceMap_()
{
// get and apply point sources if any given in the problem
@@ -1112,6 +1129,7 @@ private:
NewtonController newtonCtl_;
std::shared_ptr<VtkMultiWriter> resultWriter_;
+ std::shared_ptr<VtkOutputModule<TypeTag>> vtkOutputModule_;
std::shared_ptr<GridAdaptModel> gridAdapt_;
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
new file mode 100644
index 0000000000..0bdee3f438
--- /dev/null
+++ b/dumux/io/vtkoutputmodule.hh
@@ -0,0 +1,349 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format
+ */
+#ifndef VTK_OUTPUT_MODULE_HH
+#define VTK_OUTPUT_MODULE_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
+
+#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include <dumux/io/vtknestedfunction.hh>
+
+namespace Properties
+{
+NEW_PROP_TAG(VtkAddVelocity);
+NEW_PROP_TAG(VtkAddProcessRank);
+NEW_PROP_TAG(VtkAddPermeability);
+}
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup InputOutput
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format
+ *
+ * Handles the output of scalar and vector fields to VTK formatted file for multiple
+ * variables and timesteps. Certain predefined fields can be registered on problem / model
+ * initialization and/or be turned on/off using the designated properties. Additionally
+ * non-standardized scalar and vector fields can be added to the writer manually.
+ */
+template<typename TypeTag>
+class VtkOutputModule
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
+ using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox);
+
+ struct PriVarScalarDataInfo { unsigned int pvIdx; std::string name; };
+ struct PriVarVectorDataInfo { std::vector<unsigned int> pvIdx; std::string name; };
+ struct SecondVarScalarDataInfo { std::function<Scalar(const VolumeVariables&)> get; std::string name; };
+
+public:
+
+ VtkOutputModule(const Problem& problem,
+ Dune::VTK::DataMode dm = Dune::VTK::conforming)
+ : problem_(problem),
+ writer_(std::make_shared<Dune::VTKWriter<GridView>>(problem.gridView(), dm)),
+ sequenceWriter_(writer_, problem.name())
+ {}
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ //! Methods to conveniently add primary and secondary variables upon problem initialization
+ //! Do not call these methods after initialization
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ //! Output a scalar primary variable
+ //! \param name The name of the vtk field
+ //! \param pvIdx The index in the primary variables vector
+ void addPrimaryVariable(const std::string& name, unsigned int pvIdx)
+ {
+ priVarScalarDataInfo_.push_back(PriVarScalarDataInfo{pvIdx, name});
+ }
+
+ //! Output a vector primary variable
+ //! \param name The name of the vtk field
+ //! \param pvIndices A vector of indices in the primary variables vector to group for vector visualization
+ void addPrimaryVariable(const std::string& name, std::vector<unsigned int> pvIndices)
+ {
+ assert(pvIndices.size() < 4 && "Vtk doesn't support vector dimensions greater than 3!");
+ priVarVectorDataInfo_.push_back(PriVarVectorDataInfo{pvIndices, name});
+ }
+
+ //! Output a secondary scalar variable
+ //! \param name The name of the vtk field
+ //! \param f A function taking a VolumeVariables object and returning the desired scalar
+ void addSecondaryVariable(const std::string& name, std::function<Scalar(const VolumeVariables&)>&& f)
+ {
+ secondVarScalarDataInfo_.push_back(SecondVarScalarDataInfo{f, name});
+ }
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ //! Methods to add addtional (non-standardized) scalar or vector fields to the output queue
+ //! Call these methods on the output module obtained as argument of the addVtkOutputFields
+ //! method that can be overloaded in the problem implementation
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ //! Output a scalar field
+ //! \param name The name of the vtk field
+ //! \param codim The codimension of the entities the vtk field lives on (options: 0 := elements, dim := vertices)
+ //! \returns A reference to the resized scalar field to be filled with the actual data
+ std::vector<Scalar>& createScalarField(const std::string& name, int codim)
+ {
+ scalarFields_.emplace_back(std::make_pair(std::vector<Scalar>(problem_.gridView().size(codim)), name));
+ return scalarFields_.back().first;
+ }
+
+ //! Output a vector field
+ //! \param name The name of the vtk field
+ //! \param codim The codimension of the entities the vtk field lives on (options: 0 := elements, dim := vertices)
+ //! \returns A reference to the resized vector field to be filled with the actual data
+ std::vector<GlobalPosition>& createVectorField(const std::string& name, int codim)
+ {
+ vectorFields_.emplace_back(std::make_pair(std::vector<GlobalPosition>(problem_.gridView().size(codim)), name));
+ return vectorFields_.back().first;
+ }
+
+ //! Write the data for this timestep to file in four steps
+ //! (1) Allow user to register additional (non-standardized) fields
+ //! (2) We assemble all registered variable fields
+ //! (3) We register them with the vtk writer
+ //! (4) The writer writes the output for us
+ //! (5) Clear the writer for the next time step
+ void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
+ {
+
+ //////////////////////////////////////////////////////////////
+ //! (1) Register addtional (non-standardized) data fields with the vtk writer
+ //! Using the add scalar field or vector field methods
+ //////////////////////////////////////////////////////////////
+ problem_.addVtkOutputFields(*this);
+
+ //! Abort if no data was registered
+ //! \todo This is not necessary anymore once the old style multiwriter is removed
+ if (priVarScalarDataInfo_.empty()
+ && priVarVectorDataInfo_.empty()
+ && secondVarScalarDataInfo_.empty()
+ && scalarFields_.empty()
+ && vectorFields_.empty())
+ return;
+
+ //////////////////////////////////////////////////////////////
+ //! (2) Assemble all variable fields with registered info
+ //////////////////////////////////////////////////////////////
+ auto numCells = problem_.gridView().size(0);
+
+ // get fields for all primary variables
+ std::vector<std::vector<Scalar>> priVarScalarData(priVarScalarDataInfo_.size());
+ for (auto&& p : priVarScalarData)
+ p.resize(numCells);
+
+ std::vector<std::vector<Scalar>> priVarVectorData(priVarVectorDataInfo_.size());
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ priVarVectorData[i].resize(numCells*priVarVectorDataInfo_[i].pvIdx.size());
+
+ // get fields for all secondary variables
+ std::vector<std::vector<Scalar>> secondVarScalarData(secondVarScalarDataInfo_.size());
+ for (auto&& s : secondVarScalarData)
+ s.resize(numCells);
+
+ // instatiate the velocity output
+ ImplicitVelocityOutput<TypeTag> velocityOutput(problem_);
+ std::array<std::vector<GlobalPosition>, numPhases> velocity;
+
+ if (velocityOutput.enableOutput())
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ velocity[phaseIdx].resize(numCells);
+ }
+
+ // maybe allocate space for the process rank
+ std::vector<Scalar> rank;
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddProcessRank))
+ rank.resize(numCells);
+
+ for (const auto& element : elements(problem_.gridView(), Dune::Partitions::interior))
+ {
+ auto eIdxGlobal = problem_.elementMapper().index(element);
+
+ //! primary variable data
+ for (std::size_t i = 0; i < priVarScalarDataInfo_.size(); ++i)
+ priVarScalarData[i][eIdxGlobal] = problem_.model().curSol()[eIdxGlobal][priVarScalarDataInfo_[i].pvIdx];
+
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ for (std::size_t j = 0; j < priVarVectorDataInfo_[i].pvIdx.size(); ++j)
+ priVarVectorData[i][eIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j]
+ = problem_.model().curSol()[eIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]];
+
+ //! secondary variables
+ auto fvGeometry = localView(problem_.model().globalFvGeometry());
+ auto elemVolVars = localView(problem_.model().curGlobalVolVars());
+
+ if (velocityOutput.enableOutput() || !secondVarScalarDataInfo_.empty())
+ {
+ // If velocity output is enabled we need to bind to the whole stencil
+ // otherwise element-local data is sufficient
+ if (velocityOutput.enableOutput())
+ fvGeometry.bind(element);
+ else
+ fvGeometry.bindElement(element);
+
+ // If velocity output is enabled we need to bind to the whole stencil
+ // otherwise element-local data is sufficient
+ if (velocityOutput.enableOutput())
+ elemVolVars.bind(element, fvGeometry, problem_.model().curSol());
+ else
+ elemVolVars.bindElement(element, fvGeometry, problem_.model().curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto dofIdxGlobal = scv.dofIndex();
+ const auto& volVars = elemVolVars[scv];
+
+ for (std::size_t i = 0; i < secondVarScalarDataInfo_.size(); ++i)
+ secondVarScalarData[i][dofIdxGlobal] = secondVarScalarDataInfo_[i].get(volVars);
+ }
+ }
+
+ // velocity output
+ if (velocityOutput.enableOutput())
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
+
+ //! the rank
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddProcessRank))
+ rank[eIdxGlobal] = problem_.gridView().comm().rank();
+ }
+
+ //////////////////////////////////////////////////////////////
+ //! (3) Register data fields with the vtk writer
+ //////////////////////////////////////////////////////////////
+
+ // primary variables
+ for (std::size_t i = 0; i < priVarScalarDataInfo_.size(); ++i)
+ sequenceWriter_.addCellData(priVarScalarData[i], priVarScalarDataInfo_[i].name);
+
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ sequenceWriter_.addCellData(priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size());
+
+ // secondary variables
+ for (std::size_t i = 0; i < secondVarScalarDataInfo_.size(); ++i)
+ sequenceWriter_.addCellData(secondVarScalarData[i], secondVarScalarDataInfo_[i].name);
+
+ // the velocity field
+ if (velocityOutput.enableOutput())
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ using NestedFunction = VtkNestedFunction<GridView, ElementMapper, std::vector<GlobalPosition>>;
+ sequenceWriter_.addCellData(std::make_shared<NestedFunction>("velocity_" + std::string(FluidSystem::phaseName(phaseIdx)),
+ problem_.gridView(), problem_.elementMapper(),
+ velocity[phaseIdx], 0, dimWorld));
+ }
+ }
+
+ // the process rank
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddProcessRank))
+ sequenceWriter_.addCellData(rank, "rank");
+
+ // also register additional (non-standardized) user fields
+ for (std::size_t i = 0; i < scalarFields_.size(); ++i)
+ {
+ if (scalarFields_[i].second.size() == std::size_t(problem_.gridView().size(0)))
+ sequenceWriter_.addCellData(scalarFields_[i].first, scalarFields_[i].second);
+ else if (scalarFields_[i].second.size() == std::size_t(problem_.gridView().size(dim)))
+ sequenceWriter_.addVertexData(scalarFields_[i].first, scalarFields_[i].second);
+ else
+ DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
+ }
+
+ for (std::size_t i = 0; i < vectorFields_.size(); ++i)
+ {
+ if (scalarFields_[i].second.size() == std::size_t(problem_.gridView().size(0)))
+ {
+ using NestedFunction = VtkNestedFunction<GridView, ElementMapper, std::vector<GlobalPosition>>;
+ sequenceWriter_.addCellData(std::make_shared<NestedFunction>(vectorFields_[i].second,
+ problem_.gridView(), problem_.elementMapper(),
+ vectorFields_[i].first, 0, dimWorld));
+ }
+ else if (scalarFields_[i].second.size() == std::size_t(problem_.gridView().size(dim)))
+ {
+ using NestedFunction = VtkNestedFunction<GridView, VertexMapper, std::vector<GlobalPosition>>;
+ sequenceWriter_.addVertexData(std::make_shared<NestedFunction>(vectorFields_[i].second,
+ problem_.gridView(), problem_.vertexMapper(),
+ vectorFields_[i].first, dim, dimWorld));
+ }
+ else
+ DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk vector field!");
+ }
+
+ //////////////////////////////////////////////////////////////
+ //! (4) The writer writes the output for us
+ //////////////////////////////////////////////////////////////
+ sequenceWriter_.write(time, type);
+
+ //////////////////////////////////////////////////////////////
+ //! (5) Clear all fields and writer for the next time step
+ //////////////////////////////////////////////////////////////
+ clear();
+ }
+
+ //! clear all data in the writer
+ void clear()
+ {
+ writer_->clear();
+ scalarFields_.clear();
+ vectorFields_.clear();
+ }
+
+private:
+ const Problem& problem_;
+ std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
+ Dune::VTKSequenceWriter<GridView> sequenceWriter_;
+
+ std::vector<PriVarScalarDataInfo> priVarScalarDataInfo_;
+ std::vector<PriVarVectorDataInfo> priVarVectorDataInfo_;
+ std::vector<SecondVarScalarDataInfo> secondVarScalarDataInfo_;
+
+ std::vector<std::pair<std::vector<Scalar>, std::string>> scalarFields_;
+ std::vector<std::pair<std::vector<GlobalPosition>, std::string>> vectorFields_;
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab