24#ifndef DUMUX_IO_VTK_OUTPUT_MODULE_HH
25#define DUMUX_IO_VTK_OUTPUT_MODULE_HH
31#include <dune/common/timer.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/typetraits.hh>
35#include <dune/geometry/type.hh>
36#include <dune/geometry/multilineargeometry.hh>
38#include <dune/grid/common/mcmgmapper.hh>
39#include <dune/grid/common/partitionset.hh>
40#include <dune/grid/io/file/vtk/vtkwriter.hh>
41#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
42#include <dune/grid/common/partitionset.hh>
58template<
class Gr
idGeometry>
61 using GridView =
typename GridGeometry::GridView;
62 using Field = Vtk::template Field<GridView>;
63 static constexpr int dim = GridView::dimension;
73 const std::string&
name,
75 Dune::VTK::DataMode dm = Dune::VTK::conforming,
83 const auto precisionString = getParamFromGroup<std::string>(
paramGroup,
"Vtk.Precision",
"Float32");
86#if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
87 if (precision_ != Dumux::Vtk::Precision::float32 || coordPrecision != Dumux::Vtk::Precision::float32)
88 std::cerr <<
"Warning: Specifying VTK output precision other than Float32 is only supported in Dune 2.7 and newer. "
89 <<
"Ignoring parameter and defaulting to Float32." << std::endl;
90 writer_ = std::make_shared<Dune::VTKWriter<GridView>>(
gridGeometry.gridView(), dm);
92 writer_ = std::make_shared<Dune::VTKWriter<GridView>>(
gridGeometry.gridView(), dm, coordPrecision);
94 sequenceWriter_ = std::make_unique<Dune::VTKSequenceWriter<GridView>>(writer_,
name);
101 {
return paramGroup_; }
112 template<
typename Vector>
114 const std::string&
name,
128 template<
typename Vector>
130 const std::string&
name,
135 const auto nComp = getNumberOfComponents_(v);
137 const auto numElemDofs =
gridGeometry().elementMapper().size();
138 const auto numVertexDofs =
gridGeometry().vertexMapper().size();
143 if(numElemDofs == numVertexDofs)
144 DUNE_THROW(Dune::InvalidStateException,
"Automatic deduction of FieldType failed. Please explicitly specify FieldType::element or FieldType::vertex.");
146 if(v.size() == numElemDofs)
148 else if(v.size() == numVertexDofs)
151 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
157 if(v.size() != numElemDofs)
158 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
161 if(v.size() != numVertexDofs)
162 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
167 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.elementMapper(), v,
name, nComp, 0, dm_,
precision);
169 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.vertexMapper(), v,
name, nComp, dim, dm_,
precision);
177 void write(
double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
182 if (dm_ == Dune::VTK::conforming)
183 writeConforming_(time, type);
184 else if (dm_ == Dune::VTK::nonconforming)
185 writeNonConforming_(time, type);
187 DUNE_THROW(Dune::NotImplemented,
"Output for provided VTK data mode");
193 std::cout <<
"Writing output for problem \"" << name_ <<
"\". Took " << timer.elapsed() <<
" seconds." << std::endl;
201 const std::string&
name()
const {
return name_; }
202 Dune::VTK::DataMode
dataMode()
const {
return dm_; }
203 Dumux::Vtk::Precision
precision()
const {
return precision_; }
205 Dune::VTKWriter<GridView>&
writer() {
return *writer_; }
208 const std::vector<Field>&
fields()
const {
return fields_; }
212 virtual void writeConforming_(
double time, Dune::VTK::OutputType type)
219 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
220 std::vector<int> rank;
223 if (!fields_.empty() || addProcessRank)
225 const auto numCells = gridGeometry_.gridView().size(0);
230 rank.resize(numCells);
232 for (
const auto& element : elements(gridGeometry_.gridView(), Dune::Partitions::interior))
234 const auto eIdxGlobal = gridGeometry_.elementMapper().index(element);
235 rank[eIdxGlobal] = gridGeometry_.gridView().comm().rank();
245 this->
sequenceWriter().addCellData(Field(gridGeometry_.gridView(), gridGeometry_.elementMapper(), rank,
"process rank", 1, 0).get());
248 for (
auto&& field : fields_)
250 if (field.codim() == 0)
252 else if (field.codim() == dim)
255 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
271 virtual void writeNonConforming_(
double time, Dune::VTK::OutputType type)
273 DUNE_THROW(Dune::NotImplemented,
"Non-conforming VTK output");
277 template<
class Vector>
278 std::size_t getNumberOfComponents_(
const Vector& v)
280 if constexpr (IsIndexable<decltype(std::declval<Vector>()[0])>::value)
286 const GridGeometry& gridGeometry_;
288 const std::string paramGroup_;
289 Dune::VTK::DataMode dm_;
291 Dumux::Vtk::Precision precision_;
293 std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
294 std::unique_ptr<Dune::VTKSequenceWriter<GridView>> sequenceWriter_;
296 std::vector<Field> fields_;
311template<
class Gr
idVariables,
class SolutionVector>
315 using GridGeometry =
typename GridVariables::GridGeometry;
317 using VV =
typename GridVariables::VolumeVariables;
318 using Scalar =
typename GridVariables::Scalar;
320 using GridView =
typename GridGeometry::GridView;
323 dim = GridView::dimension,
324 dimWorld = GridView::dimensionworld
327 using Element =
typename GridView::template Codim<0>::Entity;
328 using VolVarsVector = Dune::FieldVector<Scalar, dimWorld>;
331 static constexpr int dofCodim = isBox ? dim : 0;
333 struct VolVarScalarDataInfo { std::function<Scalar(
const VV&)> get; std::string
name; Dumux::Vtk::Precision precision_; };
334 struct VolVarVectorDataInfo { std::function<VolVarsVector(
const VV&)> get; std::string
name; Dumux::Vtk::Precision precision_; };
335 using Field = Vtk::template Field<GridView>;
344 const SolutionVector&
sol,
345 const std::string&
name,
347 Dune::VTK::DataMode dm = Dune::VTK::conforming,
373 const std::string&
name)
375 volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f,
name, this->
precision()});
382 template<
class VVV = VolVarsVector,
typename std::enable_if_t<(VVV::dimension > 1),
int> = 0>
384 const std::string&
name)
386 volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f,
name, this->
precision()});
391 const auto&
problem()
const {
return gridVariables_.curGridVolVars().problem(); }
393 const GridGeometry&
gridGeometry()
const {
return gridVariables_.gridGeometry(); }
394 const SolutionVector&
sol()
const {
return sol_; }
405 void writeConforming_(
double time, Dune::VTK::OutputType type)
override
407 const Dune::VTK::DataMode dm = Dune::VTK::conforming;
414 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
417 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
418 std::vector<double> rank;
421 std::vector<std::vector<Scalar>> volVarScalarData;
422 std::vector<std::vector<VolVarsVector>> volVarVectorData;
425 if (!volVarScalarDataInfo_.empty()
426 || !volVarVectorDataInfo_.empty()
427 || !this->fields().empty()
428 || velocityOutput_->enableOutput()
431 const auto numCells =
gridGeometry().gridView().size(0);
432 const auto numDofs = numDofs_();
435 if (!volVarScalarDataInfo_.empty())
436 volVarScalarData.resize(volVarScalarDataInfo_.size(), std::vector<Scalar>(numDofs));
437 if (!volVarVectorDataInfo_.empty())
438 volVarVectorData.resize(volVarVectorDataInfo_.size(), std::vector<VolVarsVector>(numDofs));
440 if (velocityOutput_->enableOutput())
442 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
444 if(isBox && dim == 1)
445 velocity[phaseIdx].resize(numCells);
447 velocity[phaseIdx].resize(numDofs);
452 if (addProcessRank) rank.resize(numCells);
454 for (
const auto& element : elements(
gridGeometry().gridView(), Dune::Partitions::interior))
456 const auto eIdxGlobal =
gridGeometry().elementMapper().index(element);
459 auto elemVolVars =
localView(gridVariables_.curGridVolVars());
463 if (velocityOutput_->enableOutput())
465 fvGeometry.bind(element);
466 elemVolVars.bind(element, fvGeometry, sol_);
470 fvGeometry.bindElement(element);
471 elemVolVars.bindElement(element, fvGeometry, sol_);
474 if (!volVarScalarDataInfo_.empty()
475 || !volVarVectorDataInfo_.empty())
477 for (
auto&& scv : scvs(fvGeometry))
479 const auto dofIdxGlobal = scv.dofIndex();
480 const auto& volVars = elemVolVars[scv];
483 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
484 volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo_[i].get(volVars);
487 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
488 volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo_[i].get(volVars);
493 if (velocityOutput_->enableOutput())
495 auto elemFluxVarsCache =
localView(gridVariables_.gridFluxVarsCache());
496 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
498 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
499 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
504 rank[eIdxGlobal] =
static_cast<double>(
gridGeometry().gridView().comm().rank());
514 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
516 volVarScalarDataInfo_[i].
name, 1, dim, dm, this->
precision()).get() );
517 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
519 volVarVectorDataInfo_[i].
name, dimWorld, dim, dm, this->
precision()).get() );
523 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
525 volVarScalarDataInfo_[i].
name, 1, 0,dm, this->
precision()).get() );
526 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
528 volVarVectorDataInfo_[i].
name, dimWorld, 0,dm, this->
precision()).get() );
532 if (velocityOutput_->enableOutput())
534 if (isBox && dim > 1)
536 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
538 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
539 dimWorld, dim, dm, this->precision()).get() );
544 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
546 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
547 dimWorld, 0, dm, this->precision()).get() );
556 for (
auto&& field : this->
fields())
558 if (field.codim() == 0)
560 else if (field.codim() == dim)
563 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
579 void writeNonConforming_(
double time, Dune::VTK::OutputType type)
override
581 const Dune::VTK::DataMode dm = Dune::VTK::nonconforming;
584 DUNE_THROW(Dune::InvalidStateException,
"Non-conforming output makes no sense for cell-centered schemes!");
591 bool enableVelocityOutput = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddVelocity");
592 if (enableVelocityOutput ==
true && !velocityOutput_->enableOutput())
593 std::cerr <<
"Warning! Velocity output was enabled in the input file"
594 <<
" but no velocity output policy was set for the VTK output module:"
595 <<
" There will be no velocity output."
596 <<
" Use the addVelocityOutput member function of the VTK output module." << std::endl;
598 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
601 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
602 std::vector<double> rank;
605 using ScalarDataContainer = std::vector< std::vector<Scalar> >;
606 using VectorDataContainer = std::vector< std::vector<VolVarsVector> >;
607 std::vector< ScalarDataContainer > volVarScalarData;
608 std::vector< VectorDataContainer > volVarVectorData;
611 if (!volVarScalarDataInfo_.empty()
612 || !volVarVectorDataInfo_.empty()
613 || !this->fields().empty()
614 || velocityOutput_->enableOutput()
617 const auto numCells =
gridGeometry().gridView().size(0);
618 const auto numDofs = numDofs_();
621 if (!volVarScalarDataInfo_.empty())
622 volVarScalarData.resize(volVarScalarDataInfo_.size(), ScalarDataContainer(numCells));
623 if (!volVarVectorDataInfo_.empty())
624 volVarVectorData.resize(volVarVectorDataInfo_.size(), VectorDataContainer(numCells));
626 if (velocityOutput_->enableOutput())
628 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
630 if(isBox && dim == 1)
631 velocity[phaseIdx].resize(numCells);
633 velocity[phaseIdx].resize(numDofs);
638 if (addProcessRank) rank.resize(numCells);
640 for (
const auto& element : elements(
gridGeometry().gridView(), Dune::Partitions::interior))
642 const auto eIdxGlobal =
gridGeometry().elementMapper().index(element);
643 const auto numCorners = element.subEntities(dim);
646 auto elemVolVars =
localView(gridVariables_.curGridVolVars());
649 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
650 volVarScalarData[i][eIdxGlobal].resize(numCorners);
651 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
652 volVarVectorData[i][eIdxGlobal].resize(numCorners);
656 if (velocityOutput_->enableOutput())
658 fvGeometry.bind(element);
659 elemVolVars.bind(element, fvGeometry, sol_);
663 fvGeometry.bindElement(element);
664 elemVolVars.bindElement(element, fvGeometry, sol_);
667 if (!volVarScalarDataInfo_.empty()
668 || !volVarVectorDataInfo_.empty())
670 for (
auto&& scv : scvs(fvGeometry))
672 const auto& volVars = elemVolVars[scv];
675 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
676 volVarScalarData[i][scv.elementIndex()][scv.localDofIndex()] = volVarScalarDataInfo_[i].get(volVars);
679 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
680 volVarVectorData[i][scv.elementIndex()][scv.localDofIndex()] = volVarVectorDataInfo_[i].get(volVars);
685 if (velocityOutput_->enableOutput())
687 auto elemFluxVarsCache =
localView(gridVariables_.gridFluxVarsCache());
688 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
690 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
691 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
696 rank[eIdxGlobal] =
static_cast<double>(
gridGeometry().gridView().comm().rank());
704 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
706 volVarScalarDataInfo_[i].
name, 1, dim, dm, this->
precision()).get() );
708 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
710 volVarVectorDataInfo_[i].
name, dimWorld, dim, dm, this->
precision()).get() );
713 if (velocityOutput_->enableOutput())
717 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
719 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
720 dimWorld, dim, dm, this->precision()).get() );
724 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
726 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
727 dimWorld, 0,dm, this->precision()).get());
735 for (
auto&& field : this->
fields())
737 if (field.codim() == 0)
739 else if (field.codim() == dim)
742 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
758 std::size_t numDofs_()
const {
return dofCodim == dim ?
gridGeometry().vertexMapper().size() :
gridGeometry().elementMapper().size(); }
760 const GridVariables& gridVariables_;
761 const SolutionVector& sol_;
763 std::vector<VolVarScalarDataInfo> volVarScalarDataInfo_;
764 std::vector<VolVarVectorDataInfo> volVarVectorDataInfo_;
766 std::shared_ptr<VelocityOutput> velocityOutput_;
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Dune style VTK functions.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Precision stringToPrecision(std::string_view precisionName)
Maps a string (e.g. from input) to a Dune precision type.
Definition: vtkprecision.hh:53
Velocity output for implicit (porous media) models.
Definition: io/velocityoutput.hh:41
std::vector< Dune::FieldVector< Scalar, dimWorld > > VelocityVector
Definition: io/velocityoutput.hh:50
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:60
Dumux::Vtk::Precision precision() const
Definition: io/vtkoutputmodule.hh:203
const std::string & paramGroup() const
the parameter group for getting parameter from the parameter tree
Definition: io/vtkoutputmodule.hh:100
void write(double time, Dune::VTK::OutputType type=Dune::VTK::ascii)
Definition: io/vtkoutputmodule.hh:177
VtkOutputModuleBase(const GridGeometry &gridGeometry, const std::string &name, const std::string ¶mGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:72
FieldType
export field type
Definition: io/vtkoutputmodule.hh:68
Dune::VTK::DataMode dataMode() const
Definition: io/vtkoutputmodule.hh:202
void addField(const Vector &v, const std::string &name, Dumux::Vtk::Precision precision, FieldType fieldType=FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:129
Dune::VTKWriter< GridView > & writer()
Definition: io/vtkoutputmodule.hh:205
Dune::VTKSequenceWriter< GridView > & sequenceWriter()
Definition: io/vtkoutputmodule.hh:206
const std::string & name() const
Definition: io/vtkoutputmodule.hh:201
const std::vector< Field > & fields() const
Definition: io/vtkoutputmodule.hh:208
virtual ~VtkOutputModuleBase()=default
bool verbose() const
Definition: io/vtkoutputmodule.hh:200
void addField(const Vector &v, const std::string &name, FieldType fieldType=FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:113
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:198
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:313
const auto & problem() const
Definition: io/vtkoutputmodule.hh:391
const GridVariables & gridVariables() const
Definition: io/vtkoutputmodule.hh:392
void addVolumeVariable(std::function< VolVarsVector(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:383
VV VolumeVariables
export type of the volume variables for the outputfields
Definition: io/vtkoutputmodule.hh:341
void addVolumeVariable(std::function< Scalar(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:372
void addVelocityOutput(std::shared_ptr< VelocityOutputType > velocityOutput)
Add a velocity output policy.
Definition: io/vtkoutputmodule.hh:366
const VelocityOutput & velocityOutput() const
Definition: io/vtkoutputmodule.hh:400
const SolutionVector & sol() const
Definition: io/vtkoutputmodule.hh:394
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:393
VtkOutputModule(const GridVariables &gridVariables, const SolutionVector &sol, const std::string &name, const std::string ¶mGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:343
const std::vector< VolVarScalarDataInfo > & volVarScalarDataInfo() const
Definition: io/vtkoutputmodule.hh:396
const std::vector< VolVarVectorDataInfo > & volVarVectorDataInfo() const
Definition: io/vtkoutputmodule.hh:397
Velocity output for porous media models.