3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
io/vtkoutputmodule.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_IO_VTK_OUTPUT_MODULE_HH
25#define DUMUX_IO_VTK_OUTPUT_MODULE_HH
26
27#include <functional>
28#include <memory>
29#include <string>
30
31#include <dune/common/timer.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/typetraits.hh>
34
35#include <dune/geometry/type.hh>
36#include <dune/geometry/multilineargeometry.hh>
37
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>
43
44#include <dumux/common/parameters.hh>
45#include <dumux/io/format.hh>
46#include <dumux/discretization/method.hh>
47
48#include <dumux/io/vtk/function.hh>
49#include <dumux/io/vtk/fieldtype.hh>
50#include "velocityoutput.hh"
51
52namespace Dumux {
53
59template<class GridGeometry>
60class VtkOutputModuleBase
61{
62 using GridView = typename GridGeometry::GridView;
63 using Field = Vtk::template Field<GridView>;
64 static constexpr int dim = GridView::dimension;
65
66public:
68 struct [[deprecated("use Vtk::FieldType instead")]] FieldType
69 {
70 static constexpr auto element = Vtk::FieldType::element;
71 static constexpr auto vertex = Vtk::FieldType::vertex;
72 static constexpr auto automatic = Vtk::FieldType::automatic;
73 };
74
75 VtkOutputModuleBase(const GridGeometry& gridGeometry,
76 const std::string& name,
77 const std::string& paramGroup = "",
78 Dune::VTK::DataMode dm = Dune::VTK::conforming,
79 bool verbose = true)
80 : gridGeometry_(gridGeometry)
81 , name_(name)
82 , paramGroup_(paramGroup)
83 , dm_(dm)
84 , verbose_(gridGeometry.gridView().comm().rank() == 0 && verbose)
85 {
86 const auto precisionString = getParamFromGroup<std::string>(paramGroup, "Vtk.Precision", "Float32");
87 precision_ = Dumux::Vtk::stringToPrecision(precisionString);
88 const auto coordPrecision = Dumux::Vtk::stringToPrecision(getParamFromGroup<std::string>(paramGroup, "Vtk.CoordPrecision", precisionString));
89 writer_ = std::make_shared<Dune::VTKWriter<GridView>>(gridGeometry.gridView(), dm, coordPrecision);
90 sequenceWriter_ = std::make_unique<Dune::VTKSequenceWriter<GridView>>(writer_, name);
91 addProcessRank_ = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank", true);
92 }
93
94 virtual ~VtkOutputModuleBase() = default;
95
97 const std::string& paramGroup() const
98 { return paramGroup_; }
99
109 template<typename Vector>
110 void addField(const Vector& v,
111 const std::string& name,
112 Vtk::FieldType fieldType = Vtk::FieldType::automatic)
113 { addField(v, name, this->precision(), fieldType); }
114
125 template<typename Vector>
126 void addField(const Vector& v,
127 const std::string& name,
128 Dumux::Vtk::Precision precision,
129 Vtk::FieldType fieldType = Vtk::FieldType::automatic)
130 {
131 // Deduce the number of components from the given vector type
132 const auto nComp = getNumberOfComponents_(v);
133
134 const auto numElemDofs = gridGeometry().elementMapper().size();
135 const auto numVertexDofs = gridGeometry().vertexMapper().size();
136
137 // Automatically deduce the field type ...
138 if(fieldType == Vtk::FieldType::automatic)
139 {
140 if(numElemDofs == numVertexDofs)
141 DUNE_THROW(Dune::InvalidStateException, "Automatic deduction of FieldType failed. Please explicitly specify FieldType::element or FieldType::vertex.");
142
143 if(v.size() == numElemDofs)
144 fieldType = Vtk::FieldType::element;
145 else if(v.size() == numVertexDofs)
146 fieldType = Vtk::FieldType::vertex;
147 else
148 DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
149 }
150 // ... or check if the user-specified type matches the size of v
151 else
152 {
153 if(fieldType == Vtk::FieldType::element)
154 if(v.size() != numElemDofs)
155 DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
156
157 if(fieldType == Vtk::FieldType::vertex)
158 if(v.size() != numVertexDofs)
159 DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
160 }
161
162 // add the appropriate field
163 if (fieldType == Vtk::FieldType::element)
164 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.elementMapper(), v, name, nComp, 0, dm_, precision);
165 else
166 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.vertexMapper(), v, name, nComp, dim, dm_, precision);
167 }
168
174 void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
175 {
176 Dune::Timer timer;
177
178 // write to file depending on data mode
179 if (dm_ == Dune::VTK::conforming)
180 writeConforming_(time, type);
181 else if (dm_ == Dune::VTK::nonconforming)
182 writeNonConforming_(time, type);
183 else
184 DUNE_THROW(Dune::NotImplemented, "Output for provided VTK data mode");
185
187 timer.stop();
188 if (verbose_)
189 std::cout << Fmt::format("Writing output for problem \"{}\". Took {:.2g} seconds.\n", name_, timer.elapsed());
190 }
191
192protected:
193 const GridGeometry& gridGeometry() const { return gridGeometry_; }
194
195 bool verbose() const { return verbose_; }
196 const std::string& name() const { return name_; }
197 Dune::VTK::DataMode dataMode() const { return dm_; }
198 Dumux::Vtk::Precision precision() const { return precision_; }
199
200 Dune::VTKWriter<GridView>& writer() { return *writer_; }
201 Dune::VTKSequenceWriter<GridView>& sequenceWriter() { return *sequenceWriter_; }
202
203 const std::vector<Field>& fields() const { return fields_; }
204
205private:
207 virtual void writeConforming_(double time, Dune::VTK::OutputType type)
208 {
212
213 // process rank
214 std::vector<int> rank;
215
217 if (!fields_.empty() || addProcessRank_)
218 {
219 const auto numCells = gridGeometry_.gridView().size(0);
220
221 // maybe allocate space for the process rank
222 if (addProcessRank_)
223 {
224 rank.resize(numCells);
225
226 for (const auto& element : elements(gridGeometry_.gridView(), Dune::Partitions::interior))
227 {
228 const auto eIdxGlobal = gridGeometry_.elementMapper().index(element);
229 rank[eIdxGlobal] = gridGeometry_.gridView().comm().rank();
230 }
231 }
232
236
237 // the process rank
238 if (addProcessRank_)
239 this->sequenceWriter().addCellData(Field(gridGeometry_.gridView(), gridGeometry_.elementMapper(), rank, "process rank", 1, 0).get());
240
241 // also register additional (non-standardized) user fields if any
242 for (auto&& field : fields_)
243 {
244 if (field.codim() == 0)
245 this->sequenceWriter().addCellData(field.get());
246 else if (field.codim() == dim)
247 this->sequenceWriter().addVertexData(field.get());
248 else
249 DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
250 }
251 }
252
256 this->sequenceWriter().write(time, type);
257
261 this->writer().clear();
262 }
263
265 virtual void writeNonConforming_(double time, Dune::VTK::OutputType type)
266 {
267 DUNE_THROW(Dune::NotImplemented, "Non-conforming VTK output");
268 }
269
271 template<class Vector>
272 std::size_t getNumberOfComponents_(const Vector& v)
273 {
274 if constexpr (Dune::IsIndexable<decltype(std::declval<Vector>()[0])>::value)
275 return v[0].size();
276 else
277 return 1;
278 }
279
280 const GridGeometry& gridGeometry_;
281 std::string name_;
282 const std::string paramGroup_;
283 Dune::VTK::DataMode dm_;
284 bool verbose_;
285 Dumux::Vtk::Precision precision_;
286
287 std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
288 std::unique_ptr<Dune::VTKSequenceWriter<GridView>> sequenceWriter_;
289
290 std::vector<Field> fields_;
291
292 bool addProcessRank_ = true;
293};
294
307template<class GridVariables, class SolutionVector>
308class VtkOutputModule : public VtkOutputModuleBase<typename GridVariables::GridGeometry>
309{
310 using ParentType = VtkOutputModuleBase<typename GridVariables::GridGeometry>;
311 using GridGeometry = typename GridVariables::GridGeometry;
312
313 using VV = typename GridVariables::VolumeVariables;
314 using Scalar = typename GridVariables::Scalar;
315
316 using GridView = typename GridGeometry::GridView;
317
318 enum {
319 dim = GridView::dimension,
320 dimWorld = GridView::dimensionworld
321 };
322
323 using Element = typename GridView::template Codim<0>::Entity;
324 using VolVarsVector = Dune::FieldVector<Scalar, dimWorld>;
325
326 static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
327 static constexpr int dofCodim = isBox ? dim : 0;
328
329 struct VolVarScalarDataInfo { std::function<Scalar(const VV&)> get; std::string name; Dumux::Vtk::Precision precision_; };
330 struct VolVarVectorDataInfo { std::function<VolVarsVector(const VV&)> get; std::string name; Dumux::Vtk::Precision precision_; };
331 using Field = Vtk::template Field<GridView>;
332
333 using VelocityOutputType = Dumux::VelocityOutput<GridVariables>;
334
335public:
337 using VolumeVariables = VV;
338
339 VtkOutputModule(const GridVariables& gridVariables,
340 const SolutionVector& sol,
341 const std::string& name,
342 const std::string& paramGroup = "",
343 Dune::VTK::DataMode dm = Dune::VTK::conforming,
344 bool verbose = true)
345 : ParentType(gridVariables.gridGeometry(), name, paramGroup, dm, verbose)
346 , gridVariables_(gridVariables)
347 , sol_(sol)
348 , velocityOutput_(std::make_shared<VelocityOutputType>())
349 {
350 enableVelocityOutput_ = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddVelocity", false);
351 addProcessRank_ = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank", true);
352 }
353
358
365 void addVelocityOutput(std::shared_ptr<VelocityOutputType> velocityOutput)
366 { velocityOutput_ = velocityOutput; }
367
371 void addVolumeVariable(std::function<Scalar(const VolumeVariables&)>&& f,
372 const std::string& name)
373 {
374 volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f, name, this->precision()});
375 }
376
381 template<class VVV = VolVarsVector, typename std::enable_if_t<(VVV::dimension > 1), int> = 0>
382 void addVolumeVariable(std::function<VolVarsVector(const VolumeVariables&)>&& f,
383 const std::string& name)
384 {
385 volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f, name, this->precision()});
386 }
387
388protected:
389 // some return functions for differing implementations to use
390 const auto& problem() const { return gridVariables_.curGridVolVars().problem(); }
391 const GridVariables& gridVariables() const { return gridVariables_; }
392 const GridGeometry& gridGeometry() const { return gridVariables_.gridGeometry(); }
393 const SolutionVector& sol() const { return sol_; }
394
395 const std::vector<VolVarScalarDataInfo>& volVarScalarDataInfo() const { return volVarScalarDataInfo_; }
396 const std::vector<VolVarVectorDataInfo>& volVarVectorDataInfo() const { return volVarVectorDataInfo_; }
397
398 using VelocityOutput = VelocityOutputType;
399 const VelocityOutput& velocityOutput() const { return *velocityOutput_; }
400
401private:
402
404 void writeConforming_(double time, Dune::VTK::OutputType type) override
405 {
406 const Dune::VTK::DataMode dm = Dune::VTK::conforming;
410
411 // instantiate the velocity output
412 using VelocityVector = typename VelocityOutput::VelocityVector;
413 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
414
415 // process rank
416 std::vector<double> rank;
417
418 // volume variable data
419 std::vector<std::vector<Scalar>> volVarScalarData;
420 std::vector<std::vector<VolVarsVector>> volVarVectorData;
421
423 if (!volVarScalarDataInfo_.empty()
424 || !volVarVectorDataInfo_.empty()
425 || !this->fields().empty()
426 || velocityOutput_->enableOutput()
427 || addProcessRank_)
428 {
429 const auto numCells = gridGeometry().gridView().size(0);
430 const auto numDofs = numDofs_();
431
432 // get fields for all volume variables
433 if (!volVarScalarDataInfo_.empty())
434 volVarScalarData.resize(volVarScalarDataInfo_.size(), std::vector<Scalar>(numDofs));
435 if (!volVarVectorDataInfo_.empty())
436 volVarVectorData.resize(volVarVectorDataInfo_.size(), std::vector<VolVarsVector>(numDofs));
437
438 if (velocityOutput_->enableOutput())
439 {
440 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
441 {
442 if(isBox && dim == 1)
443 velocity[phaseIdx].resize(numCells);
444 else
445 velocity[phaseIdx].resize(numDofs);
446 }
447 }
448
449 // maybe allocate space for the process rank
450 if (addProcessRank_) rank.resize(numCells);
451
452 auto fvGeometry = localView(gridGeometry());
453 auto elemVolVars = localView(gridVariables_.curGridVolVars());
454
455 for (const auto& element : elements(gridGeometry().gridView(), Dune::Partitions::interior))
456 {
457 const auto eIdxGlobal = gridGeometry().elementMapper().index(element);
458 // If velocity output is enabled we need to bind to the whole stencil
459 // otherwise element-local data is sufficient
460 if (velocityOutput_->enableOutput())
461 {
462 fvGeometry.bind(element);
463 elemVolVars.bind(element, fvGeometry, sol_);
464 }
465 else
466 {
467 fvGeometry.bindElement(element);
468 elemVolVars.bindElement(element, fvGeometry, sol_);
469 }
470
471 if (!volVarScalarDataInfo_.empty()
472 || !volVarVectorDataInfo_.empty())
473 {
474 for (auto&& scv : scvs(fvGeometry))
475 {
476 const auto dofIdxGlobal = scv.dofIndex();
477 const auto& volVars = elemVolVars[scv];
478
479 // get the scalar-valued data
480 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
481 volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo_[i].get(volVars);
482
483 // get the vector-valued data
484 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
485 volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo_[i].get(volVars);
486 }
487 }
488
489 // velocity output
490 if (velocityOutput_->enableOutput())
491 {
492 const auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache()).bind(element, fvGeometry, elemVolVars);
493
494 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
495 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
496 }
497
499 if (addProcessRank_)
500 rank[eIdxGlobal] = static_cast<double>(gridGeometry().gridView().comm().rank());
501 }
502
506
507 // volume variables if any
508 if (isBox)
509 {
510 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
511 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarScalarData[i],
512 volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, dm, this->precision()).get() );
513 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
514 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarVectorData[i],
515 volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
516 }
517 else
518 {
519 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
520 this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
521 volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/0,dm, this->precision()).get() );
522 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
523 this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
524 volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/0,dm, this->precision()).get() );
525 }
526
527 // the velocity field
528 if (velocityOutput_->enableOutput())
529 {
530 if (isBox && dim > 1)
531 {
532 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
533 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
534 "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
535 /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
536 }
537 // cell-centered models
538 else
539 {
540 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
541 this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
542 "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
543 /*numComp*/dimWorld, /*codim*/0, dm, this->precision()).get() );
544 }
545 }
546
547 // the process rank
548 if (addProcessRank_)
549 this->sequenceWriter().addCellData(Field(gridGeometry().gridView(), gridGeometry().elementMapper(), rank, "process rank", 1, 0).get());
550
551 // also register additional (non-standardized) user fields if any
552 for (auto&& field : this->fields())
553 {
554 if (field.codim() == 0)
555 this->sequenceWriter().addCellData(field.get());
556 else if (field.codim() == dim)
557 this->sequenceWriter().addVertexData(field.get());
558 else
559 DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
560 }
561 }
562
566 this->sequenceWriter().write(time, type);
567
571 this->writer().clear();
572 }
573
575 void writeNonConforming_(double time, Dune::VTK::OutputType type) override
576 {
577 const Dune::VTK::DataMode dm = Dune::VTK::nonconforming;
578
579 if(!isBox)
580 DUNE_THROW(Dune::InvalidStateException, "Non-conforming output makes no sense for cell-centered schemes!");
581
585
586 // check the velocity output
587 if (enableVelocityOutput_ && !velocityOutput_->enableOutput())
588 std::cerr << "Warning! Velocity output was enabled in the input file"
589 << " but no velocity output policy was set for the VTK output module:"
590 << " There will be no velocity output."
591 << " Use the addVelocityOutput member function of the VTK output module." << std::endl;
592 using VelocityVector = typename VelocityOutput::VelocityVector;
593 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
594
595 // process rank
596 std::vector<double> rank;
597
598 // volume variable data (indexing: volvardata/element/localcorner)
599 using ScalarDataContainer = std::vector< std::vector<Scalar> >;
600 using VectorDataContainer = std::vector< std::vector<VolVarsVector> >;
601 std::vector< ScalarDataContainer > volVarScalarData;
602 std::vector< VectorDataContainer > volVarVectorData;
603
605 if (!volVarScalarDataInfo_.empty()
606 || !volVarVectorDataInfo_.empty()
607 || !this->fields().empty()
608 || velocityOutput_->enableOutput()
609 || addProcessRank_)
610 {
611 const auto numCells = gridGeometry().gridView().size(0);
612 const auto numDofs = numDofs_();
613
614 // get fields for all volume variables
615 if (!volVarScalarDataInfo_.empty())
616 volVarScalarData.resize(volVarScalarDataInfo_.size(), ScalarDataContainer(numCells));
617 if (!volVarVectorDataInfo_.empty())
618 volVarVectorData.resize(volVarVectorDataInfo_.size(), VectorDataContainer(numCells));
619
620 if (velocityOutput_->enableOutput())
621 {
622 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
623 {
624 if(isBox && dim == 1)
625 velocity[phaseIdx].resize(numCells);
626 else
627 velocity[phaseIdx].resize(numDofs);
628 }
629 }
630
631 // maybe allocate space for the process rank
632 if (addProcessRank_) rank.resize(numCells);
633
634 for (const auto& element : elements(gridGeometry().gridView(), Dune::Partitions::interior))
635 {
636 const auto eIdxGlobal = gridGeometry().elementMapper().index(element);
637 const auto numCorners = element.subEntities(dim);
638
639 auto fvGeometry = localView(gridGeometry());
640 auto elemVolVars = localView(gridVariables_.curGridVolVars());
641
642 // resize element-local data containers
643 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
644 volVarScalarData[i][eIdxGlobal].resize(numCorners);
645 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
646 volVarVectorData[i][eIdxGlobal].resize(numCorners);
647
648 // If velocity output is enabled we need to bind to the whole stencil
649 // otherwise element-local data is sufficient
650 if (velocityOutput_->enableOutput())
651 {
652 fvGeometry.bind(element);
653 elemVolVars.bind(element, fvGeometry, sol_);
654 }
655 else
656 {
657 fvGeometry.bindElement(element);
658 elemVolVars.bindElement(element, fvGeometry, sol_);
659 }
660
661 if (!volVarScalarDataInfo_.empty()
662 || !volVarVectorDataInfo_.empty())
663 {
664 for (auto&& scv : scvs(fvGeometry))
665 {
666 const auto& volVars = elemVolVars[scv];
667
668 // get the scalar-valued data
669 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
670 volVarScalarData[i][scv.elementIndex()][scv.localDofIndex()] = volVarScalarDataInfo_[i].get(volVars);
671
672 // get the vector-valued data
673 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
674 volVarVectorData[i][scv.elementIndex()][scv.localDofIndex()] = volVarVectorDataInfo_[i].get(volVars);
675 }
676 }
677
678 // velocity output
679 if (velocityOutput_->enableOutput())
680 {
681 const auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache()).bind(element, fvGeometry, elemVolVars);
682
683 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
684 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
685 }
686
688 if (addProcessRank_)
689 rank[eIdxGlobal] = static_cast<double>(gridGeometry().gridView().comm().rank());
690 }
691
695
696 // volume variables if any
697 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
698 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
699 volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, /*nonconforming*/dm, this->precision()).get() );
700
701 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
702 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
703 volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, /*nonconforming*/dm, this->precision()).get() );
704
705 // the velocity field
706 if (velocityOutput_->enableOutput())
707 {
708 // node-wise velocities
709 if (dim > 1)
710 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
711 this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
712 "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
713 /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
714
715 // cell-wise velocities
716 else
717 for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
718 this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
719 "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
720 /*numComp*/dimWorld, /*codim*/0,dm, this->precision()).get());
721 }
722
723 // the process rank
724 if (addProcessRank_)
725 this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), rank, "process rank", 1, 0).get() );
726
727 // also register additional (non-standardized) user fields if any
728 for (auto&& field : this->fields())
729 {
730 if (field.codim() == 0)
731 this->sequenceWriter().addCellData(field.get());
732 else if (field.codim() == dim)
733 this->sequenceWriter().addVertexData(field.get());
734 else
735 DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
736 }
737 }
738
742 this->sequenceWriter().write(time, type);
743
747 this->writer().clear();
748 }
749
751 std::size_t numDofs_() const { return dofCodim == dim ? gridGeometry().vertexMapper().size() : gridGeometry().elementMapper().size(); }
752
753 const GridVariables& gridVariables_;
754 const SolutionVector& sol_;
755
756 std::vector<VolVarScalarDataInfo> volVarScalarDataInfo_;
757 std::vector<VolVarVectorDataInfo> volVarVectorDataInfo_;
758
759 std::shared_ptr<VelocityOutput> velocityOutput_;
760 bool enableVelocityOutput_ = false;
761 bool addProcessRank_ = true;
762};
763
764} // end namespace Dumux
765
766#endif
parameters.hh
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
method.hh
The available discretization methods in Dumux.
fieldtype.hh
Vtk field types available in Dumux.
function.hh
Dune style VTK functions.
format.hh
Formatting based on the fmt-library which implements std::format of C++20.
Dumux::localView
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Dumux::Vtk::stringToPrecision
Precision stringToPrecision(std::string_view precisionName)
Maps a string (e.g. from input) to a Dune precision type.
Definition: precision.hh:39
Dumux::Vtk::FieldType
FieldType
Identifier for vtk field types.
Definition: fieldtype.hh:34
Dumux
Definition: adapt.hh:29
Dumux::DiscretizationMethods::box
constexpr Box box
Definition: method.hh:139
Dumux::PoreNetwork::Throat::numCorners
std::size_t numCorners(Shape shape)
Returns the number of corners of a given geometry.
Definition: throatproperties.hh:215
Dumux::VelocityOutput
Velocity output for implicit (porous media) models.
Definition: io/velocityoutput.hh:41
Dumux::VelocityOutput::VelocityVector
std::vector< Dune::FieldVector< Scalar, dimWorld > > VelocityVector
Definition: io/velocityoutput.hh:50
Dumux::VtkOutputModuleBase
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:61
Dumux::VtkOutputModuleBase::precision
Dumux::Vtk::Precision precision() const
Definition: io/vtkoutputmodule.hh:198
Dumux::VtkOutputModuleBase::paramGroup
const std::string & paramGroup() const
the parameter group for getting parameter from the parameter tree
Definition: io/vtkoutputmodule.hh:97
Dumux::VtkOutputModuleBase::write
void write(double time, Dune::VTK::OutputType type=Dune::VTK::ascii)
Definition: io/vtkoutputmodule.hh:174
Dumux::VtkOutputModuleBase::VtkOutputModuleBase
VtkOutputModuleBase(const GridGeometry &gridGeometry, const std::string &name, const std::string &paramGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:75
Dumux::VtkOutputModuleBase::dataMode
Dune::VTK::DataMode dataMode() const
Definition: io/vtkoutputmodule.hh:197
Dumux::VtkOutputModuleBase::writer
Dune::VTKWriter< GridView > & writer()
Definition: io/vtkoutputmodule.hh:200
Dumux::VtkOutputModuleBase::sequenceWriter
Dune::VTKSequenceWriter< GridView > & sequenceWriter()
Definition: io/vtkoutputmodule.hh:201
Dumux::VtkOutputModuleBase::name
const std::string & name() const
Definition: io/vtkoutputmodule.hh:196
Dumux::VtkOutputModuleBase::fields
const std::vector< Field > & fields() const
Definition: io/vtkoutputmodule.hh:203
Dumux::VtkOutputModuleBase::~VtkOutputModuleBase
virtual ~VtkOutputModuleBase()=default
Dumux::VtkOutputModuleBase::addField
void addField(const Vector &v, const std::string &name, Dumux::Vtk::Precision precision, Vtk::FieldType fieldType=Vtk::FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:126
Dumux::VtkOutputModuleBase::verbose
bool verbose() const
Definition: io/vtkoutputmodule.hh:195
Dumux::VtkOutputModuleBase::addField
void addField(const Vector &v, const std::string &name, Vtk::FieldType fieldType=Vtk::FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:110
Dumux::VtkOutputModuleBase::gridGeometry
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:193
Dumux::VtkOutputModuleBase::FieldType
export field type
Definition: io/vtkoutputmodule.hh:69
Dumux::VtkOutputModule
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:309
Dumux::VtkOutputModule::problem
const auto & problem() const
Definition: io/vtkoutputmodule.hh:390
Dumux::VtkOutputModule::gridVariables
const GridVariables & gridVariables() const
Definition: io/vtkoutputmodule.hh:391
Dumux::VtkOutputModule::addVolumeVariable
void addVolumeVariable(std::function< VolVarsVector(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:382
Dumux::VtkOutputModule::VolumeVariables
VV VolumeVariables
export type of the volume variables for the outputfields
Definition: io/vtkoutputmodule.hh:337
Dumux::VtkOutputModule::addVolumeVariable
void addVolumeVariable(std::function< Scalar(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:371
Dumux::VtkOutputModule::addVelocityOutput
void addVelocityOutput(std::shared_ptr< VelocityOutputType > velocityOutput)
Add a velocity output policy.
Definition: io/vtkoutputmodule.hh:365
Dumux::VtkOutputModule::velocityOutput
const VelocityOutput & velocityOutput() const
Definition: io/vtkoutputmodule.hh:399
Dumux::VtkOutputModule::sol
const SolutionVector & sol() const
Definition: io/vtkoutputmodule.hh:393
Dumux::VtkOutputModule::gridGeometry
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:392
Dumux::VtkOutputModule::VtkOutputModule
VtkOutputModule(const GridVariables &gridVariables, const SolutionVector &sol, const std::string &name, const std::string &paramGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:339
Dumux::VtkOutputModule::volVarScalarDataInfo
const std::vector< VolVarScalarDataInfo > & volVarScalarDataInfo() const
Definition: io/vtkoutputmodule.hh:395
Dumux::VtkOutputModule::volVarVectorDataInfo
const std::vector< VolVarVectorDataInfo > & volVarVectorDataInfo() const
Definition: io/vtkoutputmodule.hh:396
velocityoutput.hh
TODO: docme!
Generated by doxygen 1.9.3