diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index 02259ddea1b849e4c5cfdbb3fee30c77ac44fb0c..1bdbc63af13a09eae08fbc71d33253fe8e097d45 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -682,34 +682,14 @@ public:
template <class MultiWriter>
void addOutputVtkFields(const SolutionVector &sol,
MultiWriter &writer)
- {
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
-
- // create the required scalar fields
- unsigned numDofs = asImp_().numDofs();
-
- // global defect of the two auxiliary equations
- ScalarField* x[numEq];
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- x[eqIdx] = writer.allocateManagedBuffer(numDofs);
- }
-
- for (int dofIdxGlobal = 0; dofIdxGlobal < sol.size(); dofIdxGlobal++)
- {
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- (*x[eqIdx])[dofIdxGlobal] = sol[dofIdxGlobal][eqIdx];
- }
- }
+ {}
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- std::ostringstream oss;
- oss << "primaryVar_" << eqIdx;
- if (isBox)
- writer.attachVertexData(*x[eqIdx], oss.str());
- else
- writer.attachCellData(*x[eqIdx], oss.str());
- }
- }
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+ template<class VtkOutputModule>
+ void addVtkOutputFields(VtkOutputModule& outputModule) const
+ {}
/*!
* \brief Reference to the grid view of the spatial domain.
diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh
index 24675fb4af33ad390413a6f51d8c051e67a5a860..b10783c782be08354b1d130d7c82bf0bc155bee3 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/implicit/properties.hh b/dumux/implicit/properties.hh
index 8b819e2524363e7f7ff6d877a08d3430d286676c..cedbf0b2e5002e7a6b0b76b1423ef1a6ca8161c2 100644
--- a/dumux/implicit/properties.hh
+++ b/dumux/implicit/properties.hh
@@ -101,6 +101,9 @@ NEW_PROP_TAG(SolutionDependentHeatConduction); //!< specifies if the parameters
// vtk output
NEW_PROP_TAG(VtkAddVelocity); //!< specifies if an element velocity it reconstructed for the output
+NEW_PROP_TAG(VtkAddProcessRank); //!< specifies if the process rank should be added the output
+NEW_PROP_TAG(VtkAddPorosity); //!< specifies if the porosity should be added the output
+NEW_PROP_TAG(VtkAddPermeability); //!< specifies if the permeability should be added the output
// stencils
NEW_PROP_TAG(StencilsVector); //! The type of the global vector of stencils per element
diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh
index 3d6b1977b340ca1c75eaaa8667d7a97ce1eacb7a..fc55c4fa9348b585797f08aed479ac6f1e404c98 100644
--- a/dumux/implicit/propertydefaults.hh
+++ b/dumux/implicit/propertydefaults.hh
@@ -203,6 +203,9 @@ SET_SCALAR_PROP(ImplicitBase, ImplicitUpwindWeight, 1.0);
//! vtk output
SET_BOOL_PROP(ImplicitBase, VtkAddVelocity, false); //!< Don't reconstruct velocity per default
+SET_BOOL_PROP(ImplicitBase, VtkAddProcessRank, true); //!< Add process rank to output per default
+SET_BOOL_PROP(ImplicitBase, VtkAddPorosity, false); //!< Don't add porosity to output per default
+SET_BOOL_PROP(ImplicitBase, VtkAddPermeability, false); //!< Don't add permeability to output per default
} // namespace Properties
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
new file mode 100644
index 0000000000000000000000000000000000000000..088032125fca254dbabc8bb8caa2334cf2de949e
--- /dev/null
+++ b/dumux/io/vtkoutputmodule.hh
@@ -0,0 +1,412 @@
+// -*- 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(VtkAddPorosity);
+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_.model().addVtkOutputFields(*this);
+ 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);
+ auto numDofs = problem_.model().numDofs();
+
+ // get fields for all primary variables
+ std::vector<std::vector<Scalar>> priVarScalarData(priVarScalarDataInfo_.size());
+ for (auto&& p : priVarScalarData)
+ p.resize(numDofs);
+
+ std::vector<std::vector<Scalar>> priVarVectorData(priVarVectorDataInfo_.size());
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ priVarVectorData[i].resize(numDofs*priVarVectorDataInfo_[i].pvIdx.size());
+
+ // get fields for all secondary variables
+ std::vector<std::vector<Scalar>> secondVarScalarData(secondVarScalarDataInfo_.size());
+ for (auto&& s : secondVarScalarData)
+ s.resize(numDofs);
+
+ // 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(numDofs);
+ }
+
+ // maybe allocate space for the process rank
+ std::vector<Scalar> rank;
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddProcessRank))
+ rank.resize(numCells);
+
+ std::vector<Scalar> porosity;
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPorosity))
+ porosity.resize(numDofs);
+
+ std::vector<Scalar> permeability;
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPermeability))
+ permeability.resize(numDofs);
+
+ for (const auto& element : elements(problem_.gridView(), Dune::Partitions::interior))
+ {
+ const auto eIdxGlobal = problem_.elementMapper().index(element);
+
+ // cell-centered models
+ if(!isBox)
+ {
+ //! 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]];
+ }
+
+ auto fvGeometry = localView(problem_.model().globalFvGeometry());
+ auto elemVolVars = localView(problem_.model().curGlobalVolVars());
+
+ // 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();
+
+ // for box model do the privars here
+ if (isBox)
+ {
+ //! primary variable data
+ for (std::size_t i = 0; i < priVarScalarDataInfo_.size(); ++i)
+ priVarScalarData[i][dofIdxGlobal] = problem_.model().curSol()[dofIdxGlobal][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][dofIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j]
+ = problem_.model().curSol()[dofIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]];
+
+ }
+
+ // secondary variables
+ const auto& volVars = elemVolVars[scv];
+
+ for (std::size_t i = 0; i < secondVarScalarDataInfo_.size(); ++i)
+ secondVarScalarData[i][dofIdxGlobal] = secondVarScalarDataInfo_[i].get(volVars);
+
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPorosity))
+ porosity[dofIdxGlobal] = problem_.spatialParams().porosity(scv);
+
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPermeability))
+ permeability[dofIdxGlobal] = problem_.spatialParams().intrinsicPermeability(scv, 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)
+ addDofDataForWriter_(sequenceWriter_, priVarScalarData[i], priVarScalarDataInfo_[i].name);
+
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ addDofDataForWriter_(sequenceWriter_, priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size());
+
+ // secondary variables
+ for (std::size_t i = 0; i < secondVarScalarDataInfo_.size(); ++i)
+ addDofDataForWriter_(sequenceWriter_, secondVarScalarData[i], secondVarScalarDataInfo_[i].name);
+
+ // the velocity field
+ if (velocityOutput.enableOutput())
+ {
+ if (isBox)
+ {
+ using NestedFunction = VtkNestedFunction<GridView, VertexMapper, std::vector<GlobalPosition>>;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ sequenceWriter_.addVertexData(std::make_shared<NestedFunction>("velocity_" + std::string(FluidSystem::phaseName(phaseIdx)),
+ problem_.gridView(), problem_.vertexMapper(),
+ velocity[phaseIdx], dim, dimWorld));
+ }
+ // cell-centered models
+ else
+ {
+ using NestedFunction = VtkNestedFunction<GridView, ElementMapper, std::vector<GlobalPosition>>;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ 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, "process rank");
+
+ // the porosity
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPorosity))
+ addDofDataForWriter_(sequenceWriter_, porosity, "porosity");
+
+ // the permeability
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPermeability))
+ addDofDataForWriter_(sequenceWriter_, permeability, "permeability");
+
+ // also register additional (non-standardized) user fields
+ for (std::size_t i = 0; i < scalarFields_.size(); ++i)
+ {
+ if (scalarFields_[i].first.size() == std::size_t(problem_.gridView().size(0)))
+ sequenceWriter_.addCellData(scalarFields_[i].first, scalarFields_[i].second);
+ else if (scalarFields_[i].first.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].first.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].first.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:
+
+ template<typename Writer, typename... Args>
+ void addDofDataForWriter_(Writer& writer,
+ Args&&... args)
+ {
+ if (isBox)
+ writer.addVertexData(std::forward<Args>(args)...);
+ else
+ writer.addCellData(std::forward<Args>(args)...);
+ }
+
+ 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
diff --git a/dumux/material/fluidsystems/h2oairmesitylene.hh b/dumux/material/fluidsystems/h2oairmesitylene.hh
index 071c2d32d4ebd660813201829e0db2b0d0846d5c..bfe1a518c148914c46e0b52885df9456ba7960ba 100644
--- a/dumux/material/fluidsystems/h2oairmesitylene.hh
+++ b/dumux/material/fluidsystems/h2oairmesitylene.hh
@@ -187,7 +187,7 @@ public:
/*!
* \brief Return the human readable name of a phase (used in indices)
*/
- static const char *phaseName(int phaseIdx)
+ static std::string phaseName(int phaseIdx)
{
switch (phaseIdx) {
case wPhaseIdx: return "w";
@@ -200,7 +200,7 @@ public:
/*!
* \brief Return the human readable name of a component (used in indices)
*/
- static const char *componentName(int compIdx)
+ static std::string componentName(int compIdx)
{
switch (compIdx) {
case H2OIdx: return H2O::name();
diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh
index eda2fcf7b078bbe72abe573ed7b64aa06be1257f..28961f5159325cbcd9abb33534303f3e01e4779b 100644
--- a/dumux/porousmediumflow/1p/implicit/model.hh
+++ b/dumux/porousmediumflow/1p/implicit/model.hh
@@ -56,11 +56,15 @@ namespace Dumux
template<class TypeTag >
class OnePModel : public GET_PROP_TYPE(TypeTag, BaseModel)
{
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum { dim = GridView::dimension };
enum { dimWorld = GridView::dimensionworld };
@@ -68,67 +72,17 @@ class OnePModel : public GET_PROP_TYPE(TypeTag, BaseModel)
enum { dofCodim = isBox ? dim : 0 };
public:
- /*!
- * \brief \copybrief ImplicitModel::addOutputVtkFields
- *
- * Specialization for the OnePModel, adding the pressure and
- * the process rank to the VTK writer.
- */
- template<class MultiWriter>
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
+ void init(Problem& problem)
{
- // create the required scalar fields
- unsigned numDofs = this->numDofs();
+ ParentType::init(problem);
- auto &p = *(writer.allocateManagedBuffer(numDofs));
-
- auto& velocity = *(writer.template allocateManagedBuffer<double, dimWorld>(numDofs));
- ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
-
- if (velocityOutput.enableOutput())
- velocity = 0.0;
-
- unsigned numElements = this->gridView_().size(0);
- auto &rank = *(writer.allocateManagedBuffer(numElements));
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- auto eIdx = this->elementMapper().index(element);
- rank[eIdx] = this->gridView_().comm().rank();
-
- // get the local fv geometry
- auto fvGeometry = localView(this->globalFvGeometry());
- if (velocityOutput.enableOutput())
- fvGeometry.bind(element);
- else
- fvGeometry.bindElement(element);
-
- auto elemVolVars = localView(this->curGlobalVolVars());
- if (velocityOutput.enableOutput())
- elemVolVars.bind(element, fvGeometry, this->curSol());
- else
- elemVolVars.bindElement(element, fvGeometry, this->curSol());
-
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto& volVars = elemVolVars[scv];
- const auto dofIdxGlobal = scv.dofIndex();
- p[dofIdxGlobal] = volVars.pressure();
- }
-
- // velocity output
- velocityOutput.calculateVelocity(velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
- }
-
- if (velocityOutput.enableOutput())
- writer.attachDofData(velocity, "velocity", isBox, dim);
-
- writer.attachDofData(p, "p", isBox);
- writer.attachCellData(rank, "process rank");
+ // register standardized vtk output fields
+ auto& vtkOutputModule = problem.vtkOutputModule();
+ vtkOutputModule.addPrimaryVariable("pressure", Indices::pressureIdx);
}
};
-}
+
+} // end namespace Dumux
#include "propertydefaults.hh"
diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh
index 42a9e3a9e720bcd653865edf055ad8d36a569cba..b0571fb2dc6c1e54889dfc8ae7188179166ac509 100644
--- a/dumux/porousmediumflow/3p3c/implicit/model.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/model.hh
@@ -140,6 +140,48 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
enum { dofCodim = isBox ? dim : 0 };
public:
+
+ /*!
+ * \brief Apply the initial conditions to the model.
+ *
+ * \param problem The object representing the problem which needs to
+ * be simulated.
+ */
+ void init(Problem& problem)
+ {
+ ParentType::init(problem);
+
+ // register standardized vtk output fields
+ auto& vtkOutputModule = problem.vtkOutputModule();
+ vtkOutputModule.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(nPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("Sg", [](const VolumeVariables& v){ return v.saturation(gPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pg", [](const VolumeVariables& v){ return v.pressure(gPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("rhow", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("rhon", [](const VolumeVariables& v){ return v.density(nPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("rhog", [](const VolumeVariables& v){ return v.density(gPhaseIdx); });
+
+ for (int i = 0; i < numPhases; ++i)
+ for (int j = 0; j < numComponents; ++j)
+ vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(i),
+ [i,j](const VolumeVariables& v){ return v.moleFraction(i,j); });
+
+ vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+ }
+
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+ template<class VtkOutputModule>
+ void addVtkOutputFields(VtkOutputModule& outputModule) const
+ {
+ auto& phasePresence = outputModule.createScalarField("phase presence", dofCodim);
+ for (std::size_t i = 0; i < phasePresence.size(); ++i)
+ phasePresence[i] = priVarSwitch().phasePresence(i);
+ }
+
/*!
* \brief One Newton iteration was finished.
* \param uCurrent The solution after the current Newton iteration
@@ -254,137 +296,6 @@ public:
return switchFlag_;
}
- /*!
- * \brief Append all quantities of interest which can be derived
- * from the solution of the current time step to the VTK
- * writer.
- *
- * \param sol The solution vector
- * \param writer The writer for multi-file VTK datasets
- */
- template<class MultiWriter>
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
- {
- using ScalarField = Dune::BlockVector<Dune::FieldVector<double, 1> >;
- // typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
-
- // get the number of degrees of freedom
- unsigned numDofs = this->numDofs();
-
- // create the required scalar fields
- ScalarField *saturation[numPhases];
- ScalarField *pressure[numPhases];
- ScalarField *density[numPhases];
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
- {
- saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- }
-
- ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
- ScalarField *moleFraction[numPhases][numComponents];
- for (int i = 0; i < numPhases; ++i)
- for (int j = 0; j < numComponents; ++j)
- moleFraction[i][j] = writer.allocateManagedBuffer (numDofs);
- ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
- ScalarField *poro = writer.allocateManagedBuffer(numDofs);
- // VectorField *velocityN = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
- // VectorField *velocityW = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
- // VectorField *velocityG = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
- // ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
-
- // if (velocityOutput.enableOutput()) // check if velocity output is demanded
- // {
- // // initialize velocity fields
- // for (unsigned int i = 0; i < numDofs; ++i)
- // {
- // (*velocityN)[i] = Scalar(0);
- // (*velocityW)[i] = Scalar(0);
- // (*velocityG)[i] = Scalar(0);
- // }
- // }
-
- unsigned numElements = this->gridView_().size(0);
- ScalarField *rank = writer.allocateManagedBuffer (numElements);
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- int eIdx = this->problem_().elementMapper().index(element);
- (*rank)[eIdx] = this->gridView_().comm().rank();
-
- // make sure FVElementGeometry & vol vars are bound to the element
- auto fvGeometry = localView(this->globalFvGeometry());
- fvGeometry.bindElement(element);
-
- auto elemVolVars = localView(this->curGlobalVolVars());
- elemVolVars.bindElement(element, fvGeometry, this->curSol());
-
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto& volVars = elemVolVars[scv];
- int dofIdxGlobal = scv.dofIndex();
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
- {
- (*saturation[phaseIdx])[dofIdxGlobal] = volVars.saturation(phaseIdx);
- (*pressure[phaseIdx])[dofIdxGlobal] = volVars.pressure(phaseIdx);
- (*density[phaseIdx])[dofIdxGlobal] = volVars.density(phaseIdx);
- }
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- (*moleFraction[phaseIdx][compIdx])[dofIdxGlobal] = volVars.moleFraction(phaseIdx, compIdx);
-
- (*poro)[dofIdxGlobal] = volVars.porosity();
- (*temperature)[dofIdxGlobal] = volVars.temperature();
- (*phasePresence)[dofIdxGlobal] = priVarSwitch().phasePresence(dofIdxGlobal);
-
- }
-
- // velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
- // velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
- // velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
- }
-
- writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
- writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
- writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);
- writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
- writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
- writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
- writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
- writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
- writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
-
- for (int i = 0; i < numPhases; ++i)
- {
- for (int j = 0; j < numComponents; ++j)
- {
- std::ostringstream oss;
- oss << "x^"
- << FluidSystem::componentName(j)
- << "_"
- << FluidSystem::phaseName(i);
- writer.attachDofData(*moleFraction[i][j], oss.str().c_str(), isBox);
- }
- }
- writer.attachDofData(*poro, "porosity", isBox);
- writer.attachDofData(*temperature, "temperature", isBox);
- writer.attachDofData(*phasePresence, "phase presence", isBox);
-
- // if (velocityOutput.enableOutput()) // check if velocity output is demanded
- // {
- // writer.attachDofData(*velocityW, "velocityW", isBox, dim);
- // writer.attachDofData(*velocityN, "velocityN", isBox, dim);
- // writer.attachDofData(*velocityG, "velocityG", isBox, dim);
- // }
-
- writer.attachCellData(*rank, "process rank");
- }
-
/*!
* \brief Write the current solution to a restart file.
*
diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
index 186e3cebcd98d41cfe79ac9f484b68183fdb2c1f..08152b311e2315f286faf3a10ec848cde3c91788 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
@@ -327,40 +327,6 @@ public:
return wgPhaseOnly;
}
- /*!
- * \brief Append all quantities of interest which can be derived
- * from the solution of the current time step to the VTK
- * writer. Adjust this in case of anisotropic permeabilities.
- */
- void addOutputVtkFields()
- {
- // get the number of degrees of freedom
- unsigned numDofs = this->model().numDofs();
-
- // create the scalar field required for the permeabilities
- typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
- ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numDofs);
-
- for (const auto& element : elements(this->gridView()))
- {
- // make sure the FVElementGeometry is bound to the element
- auto fvGeometry = localView(this->model().globalFvGeometry());
- fvGeometry.bindElement(element);
-
- // make sure the FVElementGeometry is bound to the element
- auto elemVolVars = localView(this->model().curGlobalVolVars());
- elemVolVars.bindElement(element, fvGeometry, this->model().curSol());
-
- for (auto&& scv : scvs(fvGeometry))
- {
- auto dofIdxGlobal = scv.dofIndex();
- (*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(scv, elemVolVars[scv]);
- }
- }
-
- this->resultWriter().attachDofData(*Kxx, "permeability", isBox);
- }
-
private:
// internal method for the initial condition
void initial_(PrimaryVariables &values,
diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input b/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input
index 750bf8c828465dfbdbeba46bb1b6eee6fd75684a..9db6dcf902a4eac522f5dbdfc3280f0bc94af0bd 100644
--- a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input
+++ b/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input
@@ -12,3 +12,6 @@ Name = infiltrationcc
[Implicit]
NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences
+[Vtk]
+AddPermeability = true
+AddPorosity = true
diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input b/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input
index 553e58c86c6f5776a71cc057a3cd38cf9cec09e7..60cdfce4700b8c6a22ea88569cfddc5bb4dc45fa 100644
--- a/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input
+++ b/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input
@@ -12,3 +12,6 @@ Name = cc3p3c
[Implicit]
NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences
+[Vtk]
+AddPermeability = true
+AddPorosity = true
diff --git a/test/references/1ptestbox-reference.vtu b/test/references/1ptestbox-reference.vtu
index e6dc5207185f0c8fc61054941108aa1c02c0a884..f0873b7b0c297345b0bee2737a18fba5e7138aaf 100644
--- a/test/references/1ptestbox-reference.vtu
+++ b/test/references/1ptestbox-reference.vtu
@@ -2,8 +2,8 @@
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="100" NumberOfPoints="121">
- <PointData Scalars="p">
- <DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
+ <PointData Scalars="pressure">
+ <DataArray type="Float32" Name="pressure" NumberOfComponents="1" format="ascii">
200000 200000 191207 191581 200000 192887 200000 194783 200000 196179 200000 196611
200000 196179 200000 194783 200000 192887 200000 191581 200000 191207 181806 182275
184427 190784 193272 194034 193272 190784 184427 182275 181806 171621 171876 172210
diff --git a/test/references/1ptestcc-reference.vtu b/test/references/1ptestcc-reference.vtu
index 6e48b26ed74af2200a362070c1d30c182845e4d1..8721f71525a512d35282a846dcc0a93c6e7728b2 100644
--- a/test/references/1ptestcc-reference.vtu
+++ b/test/references/1ptestcc-reference.vtu
@@ -2,8 +2,8 @@
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="100" NumberOfPoints="121">
- <CellData Scalars="p">
- <DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
+ <CellData Scalars="pressure">
+ <DataArray type="Float32" Name="pressure" NumberOfComponents="1" format="ascii">
195910 196468 197745 198513 198838 198838 198513 197745 196468 195910 187172 188685
193744 195981 196839 196839 195981 193744 188685 187172 176921 177358 182672 186426
187727 187727 186426 182672 177358 176921 166234 166361 168478 170921 171947 171947