vtkoutputmodule.hh 35.3 KB
Newer Older
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
 * \ingroup InputOutput
22
23
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
 */
24
25
#ifndef DUMUX_IO_VTK_OUTPUT_MODULE_HH
#define DUMUX_IO_VTK_OUTPUT_MODULE_HH
26

27
#include <functional>
28
29
#include <memory>
#include <string>
30

31
#include <dune/common/timer.hh>
32
#include <dune/common/fvector.hh>
33
34
35
36
37
38
#include <dune/common/typetraits.hh>

#include <dune/geometry/type.hh>
#include <dune/geometry/multilineargeometry.hh>

#include <dune/grid/common/mcmgmapper.hh>
Katharina Heck's avatar
Katharina Heck committed
39
#include <dune/grid/common/partitionset.hh>
40
41
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
Timo Koch's avatar
Timo Koch committed
42
#include <dune/grid/common/partitionset.hh>
43

Timo Koch's avatar
Timo Koch committed
44
#include <dumux/common/parameters.hh>
45
#include <dumux/io/format.hh>
46
#include <dumux/discretization/method.hh>
47

48
#include "vtkfunction.hh"
49
#include "vtkfieldtype.hh"
Timo Koch's avatar
Timo Koch committed
50
#include "velocityoutput.hh"
Timo Koch's avatar
Timo Koch committed
51

52
namespace Dumux {
53

54
55
56
/*!
 * \ingroup InputOutput
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
57
 * \note This is a base class providing only rudimentary features
58
 */
59
60
template<class GridGeometry>
class VtkOutputModuleBase
61
{
62
    using GridView = typename GridGeometry::GridView;
63
    using Field = Vtk::template Field<GridView>;
64
    static constexpr int dim = GridView::dimension;
Timo Koch's avatar
Timo Koch committed
65

66
public:
67
    //! export field type
68
    struct [[deprecated("use Vtk::FieldType instead")]] FieldType
69
    {
70
71
72
        static constexpr auto element = Vtk::FieldType::element;
        static constexpr auto vertex = Vtk::FieldType::vertex;
        static constexpr auto automatic = Vtk::FieldType::automatic;
73
74
    };

75
76
77
78
79
80
    VtkOutputModuleBase(const GridGeometry& gridGeometry,
                        const std::string& name,
                        const std::string& paramGroup = "",
                        Dune::VTK::DataMode dm = Dune::VTK::conforming,
                        bool verbose = true)
    : 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
89
        const auto coordPrecision = Dumux::Vtk::stringToPrecision(getParamFromGroup<std::string>(paramGroup, "Vtk.CoordPrecision", precisionString));
        writer_ = std::make_shared<Dune::VTKWriter<GridView>>(gridGeometry.gridView(), dm, coordPrecision);
90
91
        sequenceWriter_ = std::make_unique<Dune::VTKSequenceWriter<GridView>>(writer_, name);
    }
92

93
    virtual ~VtkOutputModuleBase() = default;
94

95
96
97
98
    //! the parameter group for getting parameter from the parameter tree
    const std::string& paramGroup() const
    { return paramGroup_; }

Timo Koch's avatar
Timo Koch committed
99
    /*!
100
     * \brief Add a scalar or vector valued vtk field
Timo Koch's avatar
Timo Koch committed
101
     *
102
103
104
105
106
     * \param v The field to be added. Can be any indexable container. Its value type can be a number or itself an indexable container.
     * \param name The name of the field
     * \param fieldType The type of the field.
     *        This determines whether the values are associated with vertices or elements.
     *        By default, the method automatically deduces the correct type for the given input.
Timo Koch's avatar
Timo Koch committed
107
     */
108
109
110
    template<typename Vector>
    void addField(const Vector& v,
                  const std::string& name,
111
                  Vtk::FieldType fieldType = Vtk::FieldType::automatic)
112
    { addField(v, name, this->precision(), fieldType); }
113

114
115
116
117
118
119
120
121
122
123
    /*!
     * \brief Add a scalar or vector valued vtk field
     *
     * \param v The field to be added. Can be any indexable container. Its value type can be a number or itself an indexable container.
     * \param name The name of the field
     * \param fieldType The type of the field.
     *        This determines whether the values are associated with vertices or elements.
     *        By default, the method automatically deduces the correct type for the given input.
     * \param precision The output precision of this field (see Dune::VTK::Precision)
     */
124
    template<typename Vector>
125
126
    void addField(const Vector& v,
                  const std::string& name,
127
                  Dumux::Vtk::Precision precision,
128
                  Vtk::FieldType fieldType = Vtk::FieldType::automatic)
129
    {
130
131
        // Deduce the number of components from the given vector type
        const auto nComp = getNumberOfComponents_(v);
132

133
134
        const auto numElemDofs = gridGeometry().elementMapper().size();
        const auto numVertexDofs = gridGeometry().vertexMapper().size();
135
136

        // Automatically deduce the field type ...
137
        if(fieldType == Vtk::FieldType::automatic)
138
        {
139
            if(numElemDofs == numVertexDofs)
140
141
                DUNE_THROW(Dune::InvalidStateException, "Automatic deduction of FieldType failed. Please explicitly specify FieldType::element or FieldType::vertex.");

142
            if(v.size() == numElemDofs)
143
                fieldType = Vtk::FieldType::element;
144
            else if(v.size() == numVertexDofs)
145
                fieldType = Vtk::FieldType::vertex;
146
147
148
149
150
151
            else
                DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
        }
        // ... or check if the user-specified type matches the size of v
        else
        {
152
            if(fieldType == Vtk::FieldType::element)
153
                if(v.size() != numElemDofs)
154
155
                    DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");

156
            if(fieldType == Vtk::FieldType::vertex)
157
                if(v.size() != numVertexDofs)
158
159
160
161
                    DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
        }

        // add the appropriate field
162
        if (fieldType == Vtk::FieldType::element)
163
            fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.elementMapper(), v, name, nComp, 0, dm_, precision);
164
        else
165
            fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.vertexMapper(), v, name, nComp, dim, dm_, precision);
166
167
168
    }

    //! Write the data for this timestep to file in four steps
169
170
171
172
173
    //! (1) We assemble all registered variable fields
    //! (2) We register them with the vtk writer
    //! (3) The writer writes the output for us
    //! (4) Clear the writer for the next time step
    void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
174
    {
175
        Dune::Timer timer;
176

177
178
179
180
181
182
        // write to file depending on data mode
        if (dm_ == Dune::VTK::conforming)
            writeConforming_(time, type);
        else if (dm_ == Dune::VTK::nonconforming)
            writeNonConforming_(time, type);
        else
183
            DUNE_THROW(Dune::NotImplemented, "Output for provided VTK data mode");
184
185
186
187

        //! output
        timer.stop();
        if (verbose_)
188
            std::cout << Fmt::format("Writing output for problem \"{}\". Took {:.2g} seconds.\n", name_, timer.elapsed());
189
190
    }

191
protected:
192
    const GridGeometry& gridGeometry() const { return gridGeometry_; }
193

194
    bool verbose() const { return verbose_; }
195
    const std::string& name() const { return name_; }
196
    Dune::VTK::DataMode dataMode() const { return dm_; }
197
    Dumux::Vtk::Precision precision() const { return precision_; }
198
199

    Dune::VTKWriter<GridView>& writer() { return *writer_; }
200
    Dune::VTKSequenceWriter<GridView>& sequenceWriter() { return *sequenceWriter_; }
201

202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
    const std::vector<Field>& fields() const { return fields_; }

private:
    //! Assembles the fields and adds them to the writer (conforming output)
    virtual void writeConforming_(double time, Dune::VTK::OutputType type)
    {
        //////////////////////////////////////////////////////////////
        //! (1) Assemble all variable fields and add to writer
        //////////////////////////////////////////////////////////////

        // process rank
        static bool addProcessRank = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank");
        std::vector<int> rank;

       //! Abort if no data was registered
        if (!fields_.empty() || addProcessRank)
        {
            const auto numCells = gridGeometry_.gridView().size(0);

            // maybe allocate space for the process rank
            if (addProcessRank)
            {
                rank.resize(numCells);

                for (const auto& element : elements(gridGeometry_.gridView(), Dune::Partitions::interior))
                {
                    const auto eIdxGlobal = gridGeometry_.elementMapper().index(element);
                    rank[eIdxGlobal] = gridGeometry_.gridView().comm().rank();
                }
            }

            //////////////////////////////////////////////////////////////
            //! (2) Register data fields with the vtk writer
            //////////////////////////////////////////////////////////////

            // the process rank
            if (addProcessRank)
                this->sequenceWriter().addCellData(Field(gridGeometry_.gridView(), gridGeometry_.elementMapper(), rank, "process rank", 1, 0).get());

            // also register additional (non-standardized) user fields if any
            for (auto&& field : fields_)
            {
                if (field.codim() == 0)
                    this->sequenceWriter().addCellData(field.get());
                else if (field.codim() == dim)
                    this->sequenceWriter().addVertexData(field.get());
                else
                    DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
            }
        }

        //////////////////////////////////////////////////////////////
        //! (2) The writer writes the output for us
        //////////////////////////////////////////////////////////////
        this->sequenceWriter().write(time, type);

        //////////////////////////////////////////////////////////////
        //! (3) Clear the writer
        //////////////////////////////////////////////////////////////
        this->writer().clear();
    }

    //! Assembles the fields and adds them to the writer (nonconforming output)
    virtual void writeNonConforming_(double time, Dune::VTK::OutputType type)
    {
        DUNE_THROW(Dune::NotImplemented, "Non-conforming VTK output");
    }

    //! Deduces the number of components of the value type of a vector of values
    template<class Vector>
    std::size_t getNumberOfComponents_(const Vector& v)
    {
274
        if constexpr (Dune::IsIndexable<decltype(std::declval<Vector>()[0])>::value)
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
            return v[0].size();
        else
            return 1;
    }

    const GridGeometry& gridGeometry_;
    std::string name_;
    const std::string paramGroup_;
    Dune::VTK::DataMode dm_;
    bool verbose_;
    Dumux::Vtk::Precision precision_;

    std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
    std::unique_ptr<Dune::VTKSequenceWriter<GridView>> sequenceWriter_;

    std::vector<Field> fields_; //!< Registered scalar and vector fields
};

/*!
 * \ingroup InputOutput
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
 *
 * \tparam GridVariables The grid variables
 * \tparam SolutionVector The solution vector
 *
 * Handles the output of scalar and vector fields to VTK formatted file for multiple
 * variables and timesteps. Certain predefined fields can be registered on
 * 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<class GridVariables, class SolutionVector>
class VtkOutputModule : public VtkOutputModuleBase<typename GridVariables::GridGeometry>
{
    using ParentType = VtkOutputModuleBase<typename GridVariables::GridGeometry>;
    using GridGeometry = typename GridVariables::GridGeometry;

    using VV = typename GridVariables::VolumeVariables;
    using Scalar = typename GridVariables::Scalar;

    using GridView = typename GridGeometry::GridView;

    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

    using Element = typename GridView::template Codim<0>::Entity;
    using VolVarsVector = Dune::FieldVector<Scalar, dimWorld>;

    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; 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>;

public:
    //! export type of the volume variables for the outputfields
    using VolumeVariables = VV;

    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)
    : ParentType(gridVariables.gridGeometry(), name, paramGroup, dm, verbose)
    , gridVariables_(gridVariables)
    , sol_(sol)
    , velocityOutput_(std::make_shared<VelocityOutputType>())
    {}

    //////////////////////////////////////////////////////////////////////////////////////////////
    //! Methods to conveniently add primary and secondary variables upon initialization
    //! Do not call these methods after initialization i.e. _not_ within the time loop
    //////////////////////////////////////////////////////////////////////////////////////////////

    /*!
     * \brief Add a velocity output policy
     *
     * \param velocityOutput the output policy
     * \note the default policy does not add any velocity output
     */
    void addVelocityOutput(std::shared_ptr<VelocityOutputType> velocityOutput)
    { velocityOutput_ = velocityOutput; }

    //! 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)
    {
        volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f, name, this->precision()});
    }

    //! Add a vector-valued variable
    //! \param f A function taking a VolumeVariables object and returning the desired vector
    //! \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)
    {
        volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f, name, this->precision()});
    }

protected:
    // some return functions for differing implementations to use
    const auto& problem() const { return gridVariables_.curGridVolVars().problem(); }
    const GridVariables& gridVariables() const { return gridVariables_; }
    const GridGeometry& gridGeometry() const { return gridVariables_.gridGeometry(); }
    const SolutionVector& sol() const { return sol_; }

390
391
392
    const std::vector<VolVarScalarDataInfo>& volVarScalarDataInfo() const { return volVarScalarDataInfo_; }
    const std::vector<VolVarVectorDataInfo>& volVarVectorDataInfo() const { return volVarVectorDataInfo_; }

Timo Koch's avatar
Timo Koch committed
393
394
395
    using VelocityOutput = VelocityOutputType;
    const VelocityOutput& velocityOutput() const { return *velocityOutput_; }

396
397
398
private:

    //! Assembles the fields and adds them to the writer (conforming output)
399
    void writeConforming_(double time, Dune::VTK::OutputType type) override
400
    {
401
        const Dune::VTK::DataMode dm = Dune::VTK::conforming;
402
        //////////////////////////////////////////////////////////////
403
        //! (1) Assemble all variable fields and add to writer
404
405
        //////////////////////////////////////////////////////////////

406
        // instantiate the velocity output
Timo Koch's avatar
Timo Koch committed
407
        using VelocityVector = typename VelocityOutput::VelocityVector;
408
        std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
409

410
        // process rank
411
        static bool addProcessRank = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank");
412
        std::vector<double> rank;
413

414
415
        // volume variable data
        std::vector<std::vector<Scalar>> volVarScalarData;
416
        std::vector<std::vector<VolVarsVector>> volVarVectorData;
417

418
419
        //! Abort if no data was registered
        if (!volVarScalarDataInfo_.empty()
420
            || !volVarVectorDataInfo_.empty()
421
            || !this->fields().empty()
Timo Koch's avatar
Timo Koch committed
422
            || velocityOutput_->enableOutput()
423
424
            || addProcessRank)
        {
425
            const auto numCells = gridGeometry().gridView().size(0);
426
            const auto numDofs = numDofs_();
427

428
429
430
            // get fields for all volume variables
            if (!volVarScalarDataInfo_.empty())
                volVarScalarData.resize(volVarScalarDataInfo_.size(), std::vector<Scalar>(numDofs));
431
            if (!volVarVectorDataInfo_.empty())
432
                volVarVectorData.resize(volVarVectorDataInfo_.size(), std::vector<VolVarsVector>(numDofs));
433

Timo Koch's avatar
Timo Koch committed
434
            if (velocityOutput_->enableOutput())
435
            {
436
                for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
437
                {
438
439
440
441
442
443
                    if(isBox && dim == 1)
                        velocity[phaseIdx].resize(numCells);
                    else
                        velocity[phaseIdx].resize(numDofs);
                }
            }
444

445
446
            // maybe allocate space for the process rank
            if (addProcessRank) rank.resize(numCells);
447

448
            for (const auto& element : elements(gridGeometry().gridView(), Dune::Partitions::interior))
449
            {
450
                const auto eIdxGlobal = gridGeometry().elementMapper().index(element);
451

452
                auto fvGeometry = localView(gridGeometry());
453
                auto elemVolVars = localView(gridVariables_.curGridVolVars());
454

455
456
                // If velocity output is enabled we need to bind to the whole stencil
                // otherwise element-local data is sufficient
Timo Koch's avatar
Timo Koch committed
457
                if (velocityOutput_->enableOutput())
458
                {
459
                    fvGeometry.bind(element);
460
461
                    elemVolVars.bind(element, fvGeometry, sol_);
                }
462
                else
463
                {
464
465
                    fvGeometry.bindElement(element);
                    elemVolVars.bindElement(element, fvGeometry, sol_);
466
                }
467

468
469
                if (!volVarScalarDataInfo_.empty()
                    || !volVarVectorDataInfo_.empty())
470
471
472
473
474
475
                {
                    for (auto&& scv : scvs(fvGeometry))
                    {
                        const auto dofIdxGlobal = scv.dofIndex();
                        const auto& volVars = elemVolVars[scv];

476
                        // get the scalar-valued data
477
478
                        for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
                            volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo_[i].get(volVars);
479
480
481
482

                        // get the vector-valued data
                        for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
                            volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo_[i].get(volVars);
483
484
                    }
                }
485

486
                // velocity output
Timo Koch's avatar
Timo Koch committed
487
                if (velocityOutput_->enableOutput())
488
489
490
491
                {
                    auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache());
                    elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);

492
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
493
494
                        velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
                }
495

496
497
                //! the rank
                if (addProcessRank)
498
                    rank[eIdxGlobal] = static_cast<double>(gridGeometry().gridView().comm().rank());
499
            }
500

501
502
503
            //////////////////////////////////////////////////////////////
            //! (2) Register data fields with the vtk writer
            //////////////////////////////////////////////////////////////
504

505
            // volume variables if any
506
            if (isBox)
507
            {
508
                for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
509
510
                    this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarScalarData[i],
                                                         volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, dm, this->precision()).get() );
511
                for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
512
513
                    this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), volVarVectorData[i],
                                                         volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
514
515
516
517
            }
            else
            {
                for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
518
519
                    this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
                                                       volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/0,dm, this->precision()).get() );
520
                for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
521
522
                    this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
                                                       volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/0,dm, this->precision()).get() );
523
            }
524

525
            // the velocity field
Timo Koch's avatar
Timo Koch committed
526
            if (velocityOutput_->enableOutput())
527
            {
528
529
                if (isBox && dim > 1)
                {
530
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
531
                        this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
Timo Koch's avatar
Timo Koch committed
532
                                                             "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
533
                                                             /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
534
535
536
537
                }
                // cell-centered models
                else
                {
538
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
539
                        this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
Timo Koch's avatar
Timo Koch committed
540
                                                           "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
541
                                                           /*numComp*/dimWorld, /*codim*/0, dm, this->precision()).get() );
542
                }
543
544
            }

545
546
            // the process rank
            if (addProcessRank)
547
                this->sequenceWriter().addCellData(Field(gridGeometry().gridView(), gridGeometry().elementMapper(), rank, "process rank", 1, 0).get());
548

549
            // also register additional (non-standardized) user fields if any
550
            for (auto&& field : this->fields())
551
            {
552
                if (field.codim() == 0)
553
                    this->sequenceWriter().addCellData(field.get());
554
                else if (field.codim() == dim)
555
                    this->sequenceWriter().addVertexData(field.get());
556
557
                else
                    DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
558
559
560
561
            }
        }

        //////////////////////////////////////////////////////////////
562
        //! (2) The writer writes the output for us
563
        //////////////////////////////////////////////////////////////
564
        this->sequenceWriter().write(time, type);
565
566

        //////////////////////////////////////////////////////////////
567
        //! (3) Clear the writer
568
        //////////////////////////////////////////////////////////////
569
        this->writer().clear();
570
    }
571

572
    //! Assembles the fields and adds them to the writer (nonconforming output)
573
    void writeNonConforming_(double time, Dune::VTK::OutputType type) override
574
    {
575
576
        const Dune::VTK::DataMode dm = Dune::VTK::nonconforming;

577
578
579
580
581
582
583
        if(!isBox)
            DUNE_THROW(Dune::InvalidStateException, "Non-conforming output makes no sense for cell-centered schemes!");

        //////////////////////////////////////////////////////////////
        //! (1) Assemble all variable fields and add to writer
        //////////////////////////////////////////////////////////////

Timo Koch's avatar
Timo Koch committed
584
        // check the velocity output
585
        bool enableVelocityOutput = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddVelocity");
Timo Koch's avatar
Timo Koch committed
586
587
588
589
590
591
        if (enableVelocityOutput == true && !velocityOutput_->enableOutput())
            std::cerr << "Warning! Velocity output was enabled in the input file"
                      << " but no velocity output policy was set for the VTK output module:"
                      << " There will be no velocity output."
                      << " Use the addVelocityOutput member function of the VTK output module." << std::endl;
        using VelocityVector = typename VelocityOutput::VelocityVector;
592
        std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
593
594

        // process rank
595
        static bool addProcessRank = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank");
596
597
598
599
        std::vector<double> rank;

        // volume variable data (indexing: volvardata/element/localcorner)
        using ScalarDataContainer = std::vector< std::vector<Scalar> >;
600
        using VectorDataContainer = std::vector< std::vector<VolVarsVector> >;
601
602
603
604
605
606
        std::vector< ScalarDataContainer > volVarScalarData;
        std::vector< VectorDataContainer > volVarVectorData;

        //! Abort if no data was registered
        if (!volVarScalarDataInfo_.empty()
            || !volVarVectorDataInfo_.empty()
607
            || !this->fields().empty()
Timo Koch's avatar
Timo Koch committed
608
            || velocityOutput_->enableOutput()
609
            || addProcessRank)
610
        {
611
            const auto numCells = gridGeometry().gridView().size(0);
612
            const auto numDofs = numDofs_();
613
614
615
616
617
618
619

            // get fields for all volume variables
            if (!volVarScalarDataInfo_.empty())
                volVarScalarData.resize(volVarScalarDataInfo_.size(), ScalarDataContainer(numCells));
            if (!volVarVectorDataInfo_.empty())
                volVarVectorData.resize(volVarVectorDataInfo_.size(), VectorDataContainer(numCells));

Timo Koch's avatar
Timo Koch committed
620
            if (velocityOutput_->enableOutput())
621
            {
622
                for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
623
624
625
626
627
628
629
630
631
632
633
                {
                    if(isBox && dim == 1)
                        velocity[phaseIdx].resize(numCells);
                    else
                        velocity[phaseIdx].resize(numDofs);
                }
            }

            // maybe allocate space for the process rank
            if (addProcessRank) rank.resize(numCells);

634
            for (const auto& element : elements(gridGeometry().gridView(), Dune::Partitions::interior))
635
            {
636
                const auto eIdxGlobal = gridGeometry().elementMapper().index(element);
637
638
                const auto numCorners = element.subEntities(dim);

639
                auto fvGeometry = localView(gridGeometry());
640
641
642
643
644
645
646
647
648
649
                auto elemVolVars = localView(gridVariables_.curGridVolVars());

                // resize element-local data containers
                for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
                    volVarScalarData[i][eIdxGlobal].resize(numCorners);
                for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
                    volVarVectorData[i][eIdxGlobal].resize(numCorners);

                // If velocity output is enabled we need to bind to the whole stencil
                // otherwise element-local data is sufficient
Timo Koch's avatar
Timo Koch committed
650
                if (velocityOutput_->enableOutput())
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
                {
                    fvGeometry.bind(element);
                    elemVolVars.bind(element, fvGeometry, sol_);
                }
                else
                {
                    fvGeometry.bindElement(element);
                    elemVolVars.bindElement(element, fvGeometry, sol_);
                }

                if (!volVarScalarDataInfo_.empty()
                    || !volVarVectorDataInfo_.empty())
                {
                    for (auto&& scv : scvs(fvGeometry))
                    {
                        const auto& volVars = elemVolVars[scv];

                        // get the scalar-valued data
                        for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
                            volVarScalarData[i][scv.elementIndex()][scv.localDofIndex()] = volVarScalarDataInfo_[i].get(volVars);

                        // get the vector-valued data
                        for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
                            volVarVectorData[i][scv.elementIndex()][scv.localDofIndex()] = volVarVectorDataInfo_[i].get(volVars);
                    }
                }

                // velocity output
Timo Koch's avatar
Timo Koch committed
679
                if (velocityOutput_->enableOutput())
680
681
682
683
                {
                    auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache());
                    elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);

684
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
685
686
                        velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
                }
687
688
689

                //! the rank
                if (addProcessRank)
690
                    rank[eIdxGlobal] = static_cast<double>(gridGeometry().gridView().comm().rank());
691
692
693
694
695
696
697
698
            }

            //////////////////////////////////////////////////////////////
            //! Register data fields with the vtk writer
            //////////////////////////////////////////////////////////////

            // volume variables if any
            for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
699
700
                this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarScalarData[i],
                                                     volVarScalarDataInfo_[i].name, /*numComp*/1, /*codim*/dim, /*nonconforming*/dm, this->precision()).get() );
701
702

            for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
703
704
                this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), volVarVectorData[i],
                                                     volVarVectorDataInfo_[i].name, /*numComp*/dimWorld, /*codim*/dim, /*nonconforming*/dm, this->precision()).get() );
705
706

            // the velocity field
Timo Koch's avatar
Timo Koch committed
707
            if (velocityOutput_->enableOutput())
708
709
710
            {
                // node-wise velocities
                if (dim > 1)
711
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
712
                        this->sequenceWriter().addVertexData( Field(gridGeometry().gridView(), gridGeometry().vertexMapper(), velocity[phaseIdx],
Timo Koch's avatar
Timo Koch committed
713
                                                             "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
714
                                                             /*numComp*/dimWorld, /*codim*/dim, dm, this->precision()).get() );
715
716
717

                // cell-wise velocities
                else
718
                    for (int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
719
                        this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), velocity[phaseIdx],
Timo Koch's avatar
Timo Koch committed
720
                                                           "velocity_" + velocityOutput_->phaseName(phaseIdx) + " (m/s)",
721
                                                           /*numComp*/dimWorld, /*codim*/0,dm, this->precision()).get());
722
723
724
725
            }

            // the process rank
            if (addProcessRank)
726
                this->sequenceWriter().addCellData( Field(gridGeometry().gridView(), gridGeometry().elementMapper(), rank, "process rank", 1, 0).get() );
727
728

            // also register additional (non-standardized) user fields if any
729
            for (auto&& field : this->fields())
730
731
            {
                if (field.codim() == 0)
732
                    this->sequenceWriter().addCellData(field.get());
733
                else if (field.codim() == dim)
734
                    this->sequenceWriter().addVertexData(field.get());
735
736
737
                else
                    DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
            }
738
        }
739
740
741
742

        //////////////////////////////////////////////////////////////
        //! (2) The writer writes the output for us
        //////////////////////////////////////////////////////////////
743
        this->sequenceWriter().write(time, type);
744
745
746
747

        //////////////////////////////////////////////////////////////
        //! (3) Clear the writer
        //////////////////////////////////////////////////////////////
748
        this->writer().clear();
749
750
    }

751
    //! return the number of dofs, we only support vertex and cell data
752
    std::size_t numDofs_() const { return dofCodim == dim ? gridGeometry().vertexMapper().size() : gridGeometry().elementMapper().size(); }
753

754
755
756
    const GridVariables& gridVariables_;
    const SolutionVector& sol_;

757
758
759
    std::vector<VolVarScalarDataInfo> volVarScalarDataInfo_; //!< Registered volume variables (scalar)
    std::vector<VolVarVectorDataInfo> volVarVectorDataInfo_; //!< Registered volume variables (vector)

Timo Koch's avatar
Timo Koch committed
760
    std::shared_ptr<VelocityOutput> velocityOutput_; //!< The velocity output policy
761
762
763
764
765
};

} // end namespace Dumux

#endif