[vtk] Introduce new output module featuring standardized output fields

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_;
 

dumux/io/vtkoutputmodule.hh

+// -*- 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

[Dennis Gläser]

This only works for cc at the moment right?

[Dennis Gläser]

I am referring to the output module...

[Timo Koch]

yeah because maybe we want to do something more fancy for box involving svcs. Then the OutputModule has to become a property.

[Dennis Gläser]

you mean for the saturation output in case of using an interface solver? Maybe for now we can write the output such that it works for box and cc and then modify it when we have the interface solver. I suppose the output should easily work for the two methods, right? Using scv.dofIndex() to get the privars etc...

[Timo Koch]

ok I'll check

[Dennis Gläser]

Then we could merge it and do the transition to the new output module step by step... The current implementation would simply fail for a model using the box method right? There is no check whether cc is used here to skip this output no!?

[Timo Koch]

The problem is not only the interface solver. But also spatial params. The current version now continues the "bug" in the case of element parameters that randomly the last scv determines the vertex value if the parameter is not defined on the box. For solution dependent spatial params the only option is to output per scv.

[Dennis Gläser]

Of course, but what happens right now in the output module when running a box model? It should be clear and well defined before merging the branch into next...

[Timo Koch]

I guess the same happens as on master. Although the iteration sequence might be different, but I think we use the same loop of elements and local vertices.

[Timo Koch]

So on master it's not clearly defined. Do you want to improve it before or after merging? Improving it, we might need to exchange some reference solutions too.

[Dennis Gläser]

I guess I wasn't clear :). To me it seemed that the implementation of the output module only works for cc.. but is it called also for box models? You did commit calls to this writer, e.g. for the 1p model right? So what happens if I want to use the 1pbox model now.. I could also have tested it, that would've been faster I guess :). I just wanted the box output the way it was implemented here or something that makes sure we don't get a segfault when running a box model...

[Timo Koch]

I guess you didn't see the latest commit?

[Dennis Gläser]

Ah great, that's what I wanted from the beginning ;)