Commit 6fb0daa2 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/output-module' into 'next'

Feature/output module

See merge request !306
parents 719545f0 5ee0dbb6
......@@ -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.
......
......@@ -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_;
......
......@@ -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
......
......@@ -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
......
// -*- 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
......@@ -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();
......
......@@ -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);