Commit 66fc7e75 authored by Simon Emmert's avatar Simon Emmert Committed by Timo Koch
Browse files

[io] introduce possibility of variable precision for vtk output

Introduces the possibility to have variable precision, e.g. Float64.
Updates the arguments of the VTKOutputModule constructor and adds the
version check for dune. Default stays Float32.
parent 8d1fdd56
......@@ -31,6 +31,7 @@
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/io/file/vtk/function.hh>
#include <dumux/io/vtkprecision.hh>
#include <dumux/common/typetraits/typetraits.hh>
namespace Dumux {
......@@ -63,9 +64,24 @@ public:
virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>&) const
{ return accessChooser_(mycomp, mapper_.index(e), IsIndexable<decltype(field_[0])>()); }
#if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
//! get output precision for the field
Dumux::Vtk::Precision precision() const
{ return precision_; }
#else
//! get output precision for the field
Dumux::Vtk::Precision precision() const override
{ return precision_; }
#endif
//! Constructor
VectorP0VTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper)
VectorP0VTKFunction(const GridView& gridView,
const Mapper& mapper,
const F& field,
const std::string& name,
int nComps,
Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper), precision_(precision)
{
if (field.size()!=(unsigned int)(mapper.size()))
DUNE_THROW(Dune::IOError, "VectorP0VTKFunction: size mismatch between field "
......@@ -91,6 +107,7 @@ private:
const std::string name_;
int nComps_;
const Mapper& mapper_;
Dumux::Vtk::Precision precision_;
};
/*!
......@@ -131,9 +148,24 @@ public:
return interpolation.global(xi);
}
#if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
//! get output precision for the field
Dumux::Vtk::Precision precision() const
{ return precision_; }
#else
//! get output precision for the field
Dumux::Vtk::Precision precision() const override
{ return precision_; }
#endif
//! Constructor
VectorP1VTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper)
VectorP1VTKFunction(const GridView& gridView,
const Mapper& mapper,
const F& field,
const std::string& name,
int nComps,
Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper), precision_(precision)
{
if (field.size()!=(unsigned int)( mapper.size() ))
DUNE_THROW(Dune::IOError, "VectorP1VTKFunction: size mismatch between field "
......@@ -158,6 +190,7 @@ private:
const std::string name_;
int nComps_;
const Mapper& mapper_;
Dumux::Vtk::Precision precision_;
};
/*!
......@@ -202,9 +235,24 @@ public:
return interpolation.global(xi);
}
#if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
//! get output precision for the field
Dumux::Vtk::Precision precision() const
{ return precision_; }
#else
//! get output precision for the field
Dumux::Vtk::Precision precision() const override
{ return precision_; }
#endif
//! Constructor
VectorP1NonConformingVTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper)
VectorP1NonConformingVTKFunction(const GridView& gridView,
const Mapper& mapper,
const F& field,
const std::string& name,
int nComps,
Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32)
: field_(field), name_(name), nComps_(nComps), mapper_(mapper), precision_(precision)
{
if (field.size()!=(unsigned int)(mapper.size()))
DUNE_THROW(Dune::IOError, "VectorP1NonConformingVTKFunction: size mismatch between field "
......@@ -232,6 +280,8 @@ private:
const std::string name_;
int nComps_;
const Mapper& mapper_;
Dumux::Vtk::Precision precision_;
};
/*!
......@@ -252,18 +302,19 @@ public:
template <typename F, class Mapper>
Field(const GridView& gridView, const Mapper& mapper, F const& f,
const std::string& name, int numComp = 1, int codim = 0,
Dune::VTK::DataMode dm = Dune::VTK::conforming)
Dune::VTK::DataMode dm = Dune::VTK::conforming,
Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32)
: codim_(codim)
{
if (codim == GridView::dimension)
{
if (dm == Dune::VTK::conforming)
field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision);
else
field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision);
}
else if (codim == 0)
field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision);
else
DUNE_THROW(Dune::NotImplemented, "Only element or vertex quantities allowed.");
}
......@@ -276,6 +327,9 @@ public:
//! return the number of components of this field
virtual int ncomps() const { return field_->ncomps(); }
//! return the precision of this field
virtual Dumux::Vtk::Precision precision() const { return field_->precision(); }
//! codimension of the entities on which the field values live
int codim() const { return codim_; }
......@@ -296,7 +350,6 @@ private:
};
} // end namespace Vtk
} // end namespace Dumux
#endif
......@@ -82,8 +82,8 @@ class VtkOutputModule
static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
static constexpr int dofCodim = isBox ? dim : 0;
struct VolVarScalarDataInfo { std::function<Scalar(const VV&)> get; std::string name; };
struct VolVarVectorDataInfo { std::function<VolVarsVector(const VV&)> get; std::string name; };
struct VolVarScalarDataInfo { std::function<Scalar(const VV&)> get; std::string name; Dumux::Vtk::Precision precision_; };
struct VolVarVectorDataInfo { std::function<VolVarsVector(const VV&)> get; std::string name; Dumux::Vtk::Precision precision_; };
using Field = Vtk::template Field<GridView>;
using VelocityOutputType = Dumux::VelocityOutput<GridVariables>;
......@@ -108,9 +108,14 @@ public:
, sol_(sol)
, name_(name)
, paramGroup_(paramGroup)
, verbose_(gridVariables.gridGeometry().gridView().comm().rank() == 0 && verbose)
, dm_(dm)
, verbose_(gridVariables.gridGeometry().gridView().comm().rank() == 0 && verbose)
#if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
[[deprecated("Use the new VtkOutputModule with variable output precision")]]
, writer_(std::make_shared<Dune::VTKWriter<GridView>>(gridVariables.gridGeometry().gridView(), dm))
#else
, writer_(std::make_shared<Dune::VTKWriter<GridView>>(gridVariables.gridGeometry().gridView(), dm, coordPrecision_))
#endif
, sequenceWriter_(writer_, name)
, velocityOutput_(std::make_shared<VelocityOutputType>())
{}
......@@ -136,9 +141,10 @@ public:
//! Output a scalar volume variable
//! \param name The name of the vtk field
//! \param f A function taking a VolumeVariables object and returning the desired scalar
void addVolumeVariable(std::function<Scalar(const VolumeVariables&)>&& f, const std::string& name)
void addVolumeVariable(std::function<Scalar(const VolumeVariables&)>&& f,
const std::string& name)
{
volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f, name});
volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f, name, precision_});
}
//! Add a vector-valued variable
......@@ -146,9 +152,10 @@ public:
//! \param name The name of the vtk field
//! \note This method is only available for dimWorld > 1. For 1-D problems, the overload for volVar methods returning a Scalar will be used.
template<class VVV = VolVarsVector, typename std::enable_if_t<(VVV::dimension > 1), int> = 0>
void addVolumeVariable(std::function<VolVarsVector(const VolumeVariables&)>&& f, const std::string& name)
void addVolumeVariable(std::function<VolVarsVector(const VolumeVariables&)>&& f,
const std::string& name)
{
volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f, name});
volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f, name, precision_});
}
//! Add a scalar or vector valued vtk field
......@@ -158,7 +165,9 @@ public:
//! This determines whether the values are associated with vertices or elements.
//! By default, the method automatically deduces the correct type for the given input.
template<typename Vector>
void addField(const Vector& v, const std::string& name, FieldType fieldType = FieldType::automatic)
void addField(const Vector& v,
const std::string& name,
FieldType fieldType = FieldType::automatic)
{
// Deduce the number of components from the given vector type
const auto nComp = getNumberOfComponents_(v);
......@@ -193,9 +202,9 @@ public:
// add the appropriate field
if (fieldType == FieldType::element)
fields_.emplace_back(gridGeometry().gridView(), gridGeometry().elementMapper(), v, name, nComp, 0);
fields_.emplace_back(gridGeometry().gridView(), gridGeometry().elementMapper(), v, name, nComp, 0, dm_, precision_);
else
fields_.emplace_back(gridGeometry().gridView(), gridGeometry().vertexMapper(), v, name, nComp, dim);
fields_.emplace_back(gridGeometry().gridView(), gridGeometry().vertexMapper(), v, name, nComp, dim, dm_, precision_);
}
//////////////////////////////////////////////////////////////////////////////////////////////
......@@ -361,19 +370,19 @@ private:
{
for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarScalarData[i],
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim).get() );
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, dm_, precision_).get() );
for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarVectorData[i],
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim).get() );
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, dm_, precision_).get() );
}
else
{
for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
sequenceWriter_.addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/0).get() );
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/0,dm_, precision_).get() );
for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
sequenceWriter_.addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/0).get() );
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/0,dm_, precision_).get() );
}
// the velocity field
......@@ -384,7 +393,7 @@ private:
for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
"velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
/*numComp*/dimWorld, /*codim*/dim).get() );
/*numComp*/dimWorld, /*codim*/dim, dm_, precision_).get() );
}
// cell-centered models
else
......@@ -392,7 +401,7 @@ private:
for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
sequenceWriter_.addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
"velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
/*numComp*/dimWorld, /*codim*/0).get() );
/*numComp*/dimWorld, /*codim*/0, dm_, precision_).get() );
}
}
......@@ -549,11 +558,11 @@ private:
// volume variables if any
for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, /*nonconforming*/dm_).get() );
volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, /*nonconforming*/dm_, precision_).get() );
for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, /*nonconforming*/dm_).get() );
volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, /*nonconforming*/dm_, precision_).get() );
// the velocity field
if (velocityOutput_->enableOutput())
......@@ -563,14 +572,14 @@ private:
for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
sequenceWriter_.addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
"velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
/*numComp*/dimWorld, /*codim*/dim).get() );
/*numComp*/dimWorld, /*codim*/dim, dm_, precision_).get() );
// cell-wise velocities
else
for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
sequenceWriter_.addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
"velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
/*numComp*/dimWorld, /*codim*/0).get() );
/*numComp*/dimWorld, /*codim*/0,dm_, precision_).get());
}
// the process rank
......@@ -615,9 +624,12 @@ private:
const SolutionVector& sol_;
std::string name_;
std::string paramGroup_;
bool verbose_;
const std::string paramGroup_;
Dune::VTK::DataMode dm_;
bool verbose_;
const std::string precisionString_ = getParamFromGroup<std::string>(paramGroup_,"Vtk.Precision", "Float32");
const Dumux::Vtk::Precision precision_ = Dumux::Vtk::stringToPrecision(precisionString_);
const Dumux::Vtk::Precision coordPrecision_ = Dumux::Vtk::stringToPrecision(getParamFromGroup<std::string>(paramGroup_,"Vtk.CoordPrecision",precisionString_));
std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
Dune::VTKSequenceWriter<GridView> sequenceWriter_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment