diff --git a/dumux/io/vtkfunction.hh b/dumux/io/vtkfunction.hh
index 60e8e789dcb31f170a2e5bb439f662dec3b4ab1e..81ebffdffe2e829e930cdbce431dca43fe6d40be 100644
--- a/dumux/io/vtkfunction.hh
+++ b/dumux/io/vtkfunction.hh
@@ -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
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index 3584f0c9c6ed6b091380924d55cf258fdd8e8864..a6cddea5b474712a67097687ed75de3063f3f254 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -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_;