gridwriter.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_IO_GRID_WRITER_HH
#define DUMUX_IO_GRID_WRITER_HH
 
#include <config.h>
 
#if DUMUX_HAVE_GRIDFORMAT
 
#include <ranges>
#include <execution>
#include <concepts>
#include <type_traits>
 
#include <gridformat/common/type_traits.hpp>
#include <gridformat/gridformat.hpp>
#include <gridformat/traits/dune.hpp>
 
#include <dune/common/exceptions.hh>
#include <dune/common/timer.hh>
 
#include <dumux/io/format.hh>
#include <dumux/io/cvfelagrangegrid.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
 
namespace Dumux::IO {
 
 
#ifndef DOXYGEN
namespace Detail {
 
template<typename Grid, typename Format, typename Comm, typename... Args>
auto makeParallelWriter(const Grid& grid, const Format& fmt, const Dune::Communication<Comm>& comm, Args&&... args)
{
    return comm.size() > 1
        ? GridFormat::Writer<Grid>{fmt, grid, static_cast<Comm>(comm), std::forward<Args>(args)...}
        : GridFormat::Writer<Grid>{fmt, grid, std::forward<Args>(args)...};
}
 
template<typename Grid, typename Format, typename... Args>
auto makeParallelWriter(const Grid& grid, const Format& fmt, const Dune::Communication<Dune::No_Comm>&, Args&&... args)
{ return GridFormat::Writer<Grid>{fmt, grid, std::forward<Args>(args)...}; }
 
template<typename Grid, typename Format, typename... Args>
auto makeWriter(const Grid& grid, const Format& fmt, Args&&... args)
{
    const auto& comm = GridFormat::Dune::Traits::GridView<Grid>::get(grid).comm();
    return makeParallelWriter(grid, fmt, comm, std::forward<Args>(args)...);
}
 
template<typename T>
concept Container = requires(const T& t) {
    { t.size() };
    { t[std::size_t{}] };
};
 
} // namespace Detail
#endif // DOXYGEN
 
 
namespace VTK { using namespace GridFormat::VTK; }
namespace Format { using namespace GridFormat::Formats; }
namespace Encoding { using namespace GridFormat::Encoding; }
namespace Compression { using namespace GridFormat::Compression; using GridFormat::none; }
namespace Precision {
    using GridFormat::float32;
    using GridFormat::float64;
 
    using GridFormat::uint64;
    using GridFormat::uint32;
    using GridFormat::uint16;
    using GridFormat::uint8;
 
    using GridFormat::int64;
    using GridFormat::int32;
    using GridFormat::int16;
    using GridFormat::int8;
} // namespace Precision
 

template<int order>
struct Order { static_assert(order > 0, "order must be > 0"); };
 
template<int o>
inline constexpr auto order = Order<o>{};

template<GridFormat::Concepts::Grid GridView, int order = 1>
class GridWriter
{
    using Grid = std::conditional_t<
        (order > 1),
        CVFELagrangeGrid<GridView, order>,
        GridView
    >;
    using Cell = GridFormat::Cell<Grid>;
    using Vertex = typename GridView::template Codim<GridView::dimension>::Entity;
    using Element = typename GridView::template Codim<0>::Entity;
    using Writer = GridFormat::Writer<Grid>;
 
 public:
    template<typename Format>
    explicit GridWriter(const Format& fmt,
                        const GridView& gridView,
                        const Order<order>& = {})
    : gridView_{gridView}
    , grid_{makeGrid_(gridView)}
    , writer_{Detail::makeWriter(grid_, fmt)}
    {
        if (gridView.comm().size() > 0 && getParam<bool>("IO.GridWriter.AddProcessRank", true))
            setCellField("process rank", [&](const Cell&) { return gridView.comm().rank(); }, Precision::uint64);
    }

    template<typename Format>
    explicit GridWriter(const Format& fmt,
                        const GridView& gridView,
                        const std::string& filename,
                        const Order<order>& = {})
    : gridView_{gridView}
    , grid_{makeGrid_(gridView)}
    , writer_{Detail::makeWriter(grid_, fmt, filename)}
    {
        if (gridView.comm().size() > 0 && getParam<bool>("IO.GridWriter.AddProcessRank", true))
            setCellField("process rank", [&](const Cell&) { return gridView.comm().rank(); }, Precision::uint64);
    }

    std::string write(const std::string& name) const
    { return writer_.write(name); }

    template<std::floating_point T>
    std::string write(T time) const
    { return writer_.write(time); }

    template<GridFormat::Concepts::CellFunction<GridView> F,
             GridFormat::Concepts::Scalar T = GridFormat::FieldScalar<std::invoke_result_t<F, Element>>>
    void setCellField(const std::string& name, F&& f, const GridFormat::Precision<T>& prec = {})
    { writer_.set_cell_field(name, std::move(f), prec); }

    template<GridFormat::Dune::Concepts::Function<GridView> F>
    void setCellField(const std::string& name, F&& f)
    { GridFormat::Dune::set_cell_function(std::forward<F>(f), writer_, name); }

    template<typename F, GridFormat::Concepts::Scalar T>
    void setCellField(const std::string& name, F&& f, const GridFormat::Precision<T>& prec)
    { GridFormat::Dune::set_cell_function(std::forward<F>(f), writer_, name, prec); }

    template<GridFormat::Concepts::PointFunction<GridView> F,
             GridFormat::Concepts::Scalar T = GridFormat::FieldScalar<std::invoke_result_t<F, Vertex>>>
    void setPointField(const std::string& name, F&& f, const GridFormat::Precision<T>& prec = {})
        requires (order == 1)
    { writer_.set_point_field(name, std::move(f), prec); }

    template<Detail::Container C,
             GridFormat::Concepts::Scalar T = GridFormat::MDRangeScalar<C>>
    void setPointField(const std::string& name, const C& values,
                       const GridFormat::Precision<T>& prec = {})
        requires (order > 1)
    {
        using Point = typename Grid::Point;
        writer_.set_point_field(name, [&values](const Point& p) -> std::ranges::range_value_t<C> {
            return values[p.dofIndex];
        }, prec);
    }

    template<GridFormat::Dune::Concepts::Function<GridView> F>
    void setPointField(const std::string& name, F&& f)
    { GridFormat::Dune::set_point_function(std::forward<F>(f), writer_, name); }

    template<GridFormat::Dune::Concepts::Function<GridView> F, GridFormat::Concepts::Scalar T>
    void setPointField(const std::string& name, F&& f, const GridFormat::Precision<T>& prec)
    { GridFormat::Dune::set_point_function(std::forward<F>(f), writer_, name, prec); }

    void clear()
    { writer_.clear(); }

    void update()
    {
        if constexpr (order > 1)
            grid_.update(gridView_);
    }
 
 private:
    Grid makeGrid_(const GridView& gv) const
    {
        if constexpr (order > 1)
            return Grid{gv};
        else
            return gv;
    }
 
    GridView gridView_;
    Grid grid_;
    Writer writer_;
};

template<typename GridVariables, typename SolutionVector>
class OutputModule : private GridWriter<typename GridVariables::GridGeometry::GridView, 1> {
    using ParentType = GridWriter<typename GridVariables::GridGeometry::GridView, 1>;
    using GridView = typename GridVariables::GridGeometry::GridView;
 
    static constexpr bool isCVFE = DiscretizationMethods::isCVFE<typename GridVariables::GridGeometry::DiscretizationMethod>;
    static constexpr int dimWorld = GridVariables::GridGeometry::GridView::dimensionworld;
    using Scalar = typename GridVariables::Scalar;
    using Vector = Dune::FieldVector<Scalar, dimWorld>;
    using Tensor = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
    using VolVar = typename GridVariables::VolumeVariables;
 
    class VolVarFieldStorage;
 
    static constexpr int defaultVerbosity_ = 1;
 
 public:
    using VolumeVariables = VolVar;
 
    static constexpr auto defaultFileFormat = IO::Format::pvd_with(
        IO::Format::vtu.with({
            .encoder = IO::Encoding::ascii,
            .compressor = IO::Compression::none,
            .data_format = VTK::DataFormat::inlined
        })
    );

    explicit OutputModule(const GridVariables& gridVariables,
                          const SolutionVector& sol,
                          const std::string& filename)
    : ParentType{defaultFileFormat, gridVariables.gridGeometry().gridView(), filename, order<1>}
    , gridVariables_{gridVariables}
    , solutionVector_{sol}
    {
        setVerbosity(defaultVerbosity_);
    }

    template<typename Format>
    explicit OutputModule(const Format& fmt,
                          const GridVariables& gridVariables,
                          const SolutionVector& sol)
    : ParentType{fmt, gridVariables.gridGeometry().gridView(), order<1>}
    , gridVariables_{gridVariables}
    , solutionVector_{sol}
    {
        setVerbosity(defaultVerbosity_);
    }

    template<typename Format>
    explicit OutputModule(const Format& fmt,
                          const GridVariables& gridVariables,
                          const SolutionVector& sol,
                          const std::string& filename)
    : ParentType{fmt, gridVariables.gridGeometry().gridView(), filename, order<1>}
    , gridVariables_{gridVariables}
    , solutionVector_{sol}
    {
        setVerbosity(defaultVerbosity_);
    }

    template<std::invocable<const VolumeVariables&> VolVarFunction>
    void addVolumeVariable(VolVarFunction&& f, const std::string& name)
    {
        using ResultType = std::invoke_result_t<std::remove_cvref_t<VolVarFunction>, const VolumeVariables&>;
        if constexpr (GridFormat::Concepts::Scalar<ResultType>)
            setVolVarField_<ResultType>(name, volVarFields_.registerScalarField(name, [_f=std::move(f)] (const auto& vv) {
                return static_cast<Scalar>(_f(vv));
            }));
        else if constexpr (GridFormat::mdrange_dimension<ResultType> == 1)
            setVolVarField_<GridFormat::MDRangeScalar<ResultType>>(name, volVarFields_.registerVectorField(name, [_f=std::move(f)] (const auto& vv) {
                return VolVarFieldStorage::toStorageVector(_f(vv));
            }));
        else if constexpr (GridFormat::mdrange_dimension<ResultType> == 2)
            setVolVarField_<GridFormat::MDRangeScalar<ResultType>>(name, volVarFields_.registerTensorField(name, [_f=std::move(f)] (const auto& vv) {
                return VolVarFieldStorage::toStorageTensor(_f(vv));
            }));
        else
        {
            static_assert(
                Dune::AlwaysFalse<VolVarFunction>::value,
                "Could not identify the given volume variable as scalar, vector or tensor."
            );
        }
    }

    template<Detail::Container C>
    void addField(const C& values, const std::string& name)
    { addField(values, name, GridFormat::Precision<GridFormat::MDRangeScalar<C>>{}); }

    template<Detail::Container C, GridFormat::Concepts::Scalar T>
    void addField(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        const bool hasCellSize = values.size() == gridVariables_.gridGeometry().elementMapper().size();
        const bool hasPointSize = values.size() == gridVariables_.gridGeometry().vertexMapper().size();
        if (hasCellSize && hasPointSize)
            DUNE_THROW(Dune::InvalidStateException, "Automatic deduction of field type failed. Please use addCellField or addPointField instead.");
        if (!hasCellSize && !hasPointSize)
            DUNE_THROW(Dune::InvalidStateException, "Automatic deduction of field type failed. Given container size does not match neither the number of points nor cells.");
 
        if (hasCellSize)
            addCellField(values, name, prec);
        else
            addPointField(values, name, prec);
    }

    template<GridFormat::Concepts::PointFunction<GridView> DofFunction>
    void addPointField(DofFunction&& f, const std::string& name)
    { this->setPointField(name, std::forward<DofFunction>(f)); }

    template<Detail::Container C>
    void addPointField(const C& values, const std::string& name)
    { addPointField(values, name, GridFormat::Precision<GridFormat::MDRangeScalar<C>>{}); }

    template<Detail::Container C, GridFormat::Concepts::Scalar T>
    void addPointField(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        if (values.size() != gridVariables_.gridGeometry().vertexMapper().size())
            DUNE_THROW(Dune::InvalidStateException, "Given container does not match the number of points in the grid");
 
        addPointField_(values, name, prec);
    }

    template<GridFormat::Concepts::CellFunction<GridView> DofFunction>
    void addCellField(DofFunction&& f, const std::string& name)
    { this->setCellField(name, std::forward<DofFunction>(f)); }

    template<Detail::Container C>
    void addCellField(const C& values, const std::string& name)
    { addCellField(values, name, GridFormat::Precision<GridFormat::MDRangeScalar<C>>{}); }

    template<Detail::Container C, GridFormat::Concepts::Scalar T>
    void addCellField(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        if (values.size() != gridVariables_.gridGeometry().elementMapper().size())
            DUNE_THROW(Dune::InvalidStateException, "Given container does not match the number of cells in the grid");
 
        addCellField_(values, name, prec);
    }

    std::string write(const std::string& name)
    {
        Dune::Timer timer;
 
        volVarFields_.updateFieldData(gridVariables_, solutionVector_);
        auto filename = ParentType::write(name);
        volVarFields_.clearFieldData();
 
        timer.stop();
        if (verbosity_ > 0)
            std::cout << Fmt::format(
                "Writing output to \"{}\". Took {:.2g} seconds.\n",
                filename, timer.elapsed()
            );
 
        return filename;
    }

    template<std::floating_point T>
    std::string write(T time)
    {
        Dune::Timer timer;
 
        volVarFields_.updateFieldData(gridVariables_, solutionVector_);
        auto filename = ParentType::write(time);
        volVarFields_.clearFieldData();
 
        timer.stop();
        if (verbosity_ > 0)
            std::cout << Fmt::format(
                "Writing output to \"{}\". Took {:.2g} seconds.\n",
                filename, timer.elapsed()
            );
 
        return filename;
     }

    void clear() {
        ParentType::clear();
        volVarFields_.clear();
    }
 
    void setVerbosity(int verbosity)
    {
        if (gridVariables_.gridGeometry().gridView().comm().rank() == 0)
            verbosity_ = verbosity;
    }
 
private:
    // reproduce behaviour of VtkOutputModule that flattens vectors of size 1
    template<Detail::Container C, GridFormat::Concepts::Scalar T>
        requires(GridFormat::has_sub_range<C> && std::ranges::range_value_t<C>::size() == 1)
    void addPointField_(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        this->setPointField(name, [&] (const auto& vertex) -> T {
            return values[gridVariables_.gridGeometry().vertexMapper().index(vertex)][0];
        }, prec);
    }
 
    template<Detail::Container C, GridFormat::Concepts::Scalar T>
    void addPointField_(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        this->setPointField(name, [&] (const auto& vertex) -> std::ranges::range_value_t<C> {
            return values[gridVariables_.gridGeometry().vertexMapper().index(vertex)];
        }, prec);
    }
 
    // reproduce behaviour of VtkOutputModule that flattens vectors of size 1
    template<Detail::Container C, GridFormat::Concepts::Scalar T>
        requires(GridFormat::has_sub_range<C> && std::ranges::range_value_t<C>::size() == 1)
    void addCellField_(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        this->setCellField(name, [&] (const auto& element) -> T {
            return values[gridVariables_.gridGeometry().elementMapper().index(element)][0];
        }, prec);
    }
 
    template<Detail::Container C, GridFormat::Concepts::Scalar T>
    void addCellField_(const C& values, const std::string& name, const GridFormat::Precision<T>& prec)
    {
        this->setCellField(name, [&] (const auto& element) -> std::ranges::range_value_t<C> {
            return values[gridVariables_.gridGeometry().elementMapper().index(element)];
        }, prec);
    }
 
    template<typename ResultType, typename Id>
    void setVolVarField_(const std::string& name, Id&& volVarFieldId)
    {
        auto dofEntityField = [&, _id=std::move(volVarFieldId)] (const auto& entity) {
            return volVarFields_.getValue(_id, gridVariables_.gridGeometry().dofMapper().index(entity));
        };
        if constexpr (isCVFE)
            this->setPointField(name, std::move(dofEntityField), GridFormat::Precision<ResultType>{});
        else
            this->setCellField(name, std::move(dofEntityField), GridFormat::Precision<ResultType>{});
    }
 
    const GridVariables& gridVariables_;
    const SolutionVector& solutionVector_;
    VolVarFieldStorage volVarFields_;
    int verbosity_ = 0;
};
 
// Class to store vol var fields; implementation detail of the OutputModule class
template<typename GridVariables, typename SolutionVector>
class OutputModule<GridVariables, SolutionVector>::VolVarFieldStorage
{
    enum class FieldType
    { scalar, vector, tensor };
 
    struct FieldInfo
    {
        std::string name;
        FieldType type;
        std::size_t index;
    };
 
    template<typename T>
    struct FieldStorage
    {
        std::vector<T> data;
        std::function<T(const VolVar&)> getter;
    };
 
public:
    template<FieldType ft>
    struct FieldId { std::size_t index; };
 
    template<std::ranges::range R>
    static constexpr auto toStorageVector(R&& in)
    {
        Vector result;
        std::ranges::copy(in, result.begin());
        return result;
    }
 
    template<GridFormat::Concepts::MDRange<2> R>
    static constexpr auto toStorageTensor(R&& in)
    {
        Tensor result;
        std::ranges::for_each(in, [&, i=0] (const auto& row) mutable {
            std::ranges::copy(row, result[i++].begin());
        });
        return result;
    }
 
    template<FieldType ft>
    const auto& getValue(const FieldId<ft>& id, std::size_t idx) const
    {
        if constexpr (ft == FieldType::scalar)
            return scalarFieldStorage_.at(id.index).data.at(idx);
        else if constexpr (ft == FieldType::vector)
            return vectorFieldStorage_.at(id.index).data.at(idx);
        else
            return tensorFieldStorage_.at(id.index).data.at(idx);
    }
 
    auto registerScalarField(std::string name, std::function<Scalar(const VolVar&)> f)
    { return register_<FieldType::scalar>(std::move(name), scalarFieldStorage_, std::move(f)); }
 
    auto registerVectorField(std::string name, std::function<Vector(const VolVar&)> f)
    { return register_<FieldType::vector>(std::move(name), vectorFieldStorage_, std::move(f)); }
 
    auto registerTensorField(std::string name, std::function<Tensor(const VolVar&)> f)
    { return register_<FieldType::tensor>(std::move(name), tensorFieldStorage_, std::move(f)); }
 
    void updateFieldData(const GridVariables& gridVars, const SolutionVector& x)
    {
        resizeFieldData_(gridVars.gridGeometry().numDofs());
        const auto range = GridFormat::cells(gridVars.gridGeometry().gridView());
        std::for_each(
#if __cpp_lib_parallel_algorithm >= 201603L
            std::execution::par_unseq,
#endif
            std::ranges::begin(range),
            std::ranges::end(range),
            [&] (const auto& element) {
                auto fvGeometry = localView(gridVars.gridGeometry()).bindElement(element);
                auto elemVolVars = localView(gridVars.curGridVolVars()).bindElement(element, fvGeometry, x);
                for (const auto& scv : scvs(fvGeometry))
                {
                    const auto& volVars = elemVolVars[scv];
                    for (auto& s : scalarFieldStorage_) { s.data.at(scv.dofIndex()) = s.getter(volVars); }
                    for (auto& s : vectorFieldStorage_) { s.data.at(scv.dofIndex()) = s.getter(volVars); }
                    for (auto& s : tensorFieldStorage_) { s.data.at(scv.dofIndex()) = s.getter(volVars); }
                }
        });
    }
 
    void clearFieldData()
    {
        for (auto& s : scalarFieldStorage_) { s.data.clear(); }
        for (auto& s : vectorFieldStorage_) { s.data.clear(); }
        for (auto& s : tensorFieldStorage_) { s.data.clear(); }
    }
 
    void clear()
    {
        fields_.clear();
        scalarFieldStorage_.clear();
        vectorFieldStorage_.clear();
        tensorFieldStorage_.clear();
    }
 
private:
    template<FieldType ft, typename T>
    auto register_(std::string&& name,
                   std::vector<FieldStorage<T>>& storage,
                   std::function<T(const VolVar&)>&& f)
    {
        if (exists_<ft>(name))
            DUNE_THROW(Dune::InvalidStateException, "Volume variables field '" << name << "' is already defined.");
 
        FieldId<ft> id{storage.size()};
        fields_.emplace_back(FieldInfo{std::move(name), ft, id.index});
        storage.push_back({{}, std::move(f)});
        return id;
    }
 
    template<FieldType ft>
    bool exists_(const std::string& name) const
    {
        return std::ranges::any_of(fields_, [&] (const FieldInfo& info) {
            return info.type == ft && info.name == name;
        });
    }
 
    void resizeFieldData_(std::size_t size)
    {
        std::ranges::for_each(scalarFieldStorage_, [&] (auto& s) { s.data.resize(size); });
        std::ranges::for_each(vectorFieldStorage_, [&] (auto& s) { s.data.resize(size); });
        std::ranges::for_each(tensorFieldStorage_, [&] (auto& s) { s.data.resize(size); });
    }
 
    std::vector<FieldInfo> fields_;
    std::vector<FieldStorage<Scalar>> scalarFieldStorage_;
    std::vector<FieldStorage<Vector>> vectorFieldStorage_;
    std::vector<FieldStorage<Tensor>> tensorFieldStorage_;
};
 
} // namespace Dumux::IO
 
#else // DUMUX_HAVE_GRIDFORMAT
 
namespace Dumux::IO {
 
template<class... Args>
class GridWriter
{
public:
    template<class... _Args>
    GridWriter(_Args&&...)
    {
        static_assert(
            false,
            "GridWriter only available when the GridFormat library is available. "
            "Use `git submodule update --init` to pull it and reconfigure the project "
            "(note: C++20 is required)."
        );
    }
};
 
template<class... Args>
class OutputModule
{
public:
    template<class... _Args>
    OutputModule(_Args&&...)
    {
        static_assert(
            false,
            "OutputModule only available when the GridFormat library is available. "
            "Use `git submodule update --init` to pull it and reconfigure the project "
            "(note: C++20 is required)."
        );
    }
};
 
} // namespace Dumux::IO
 
#endif // DUMUX_HAVE_GRIDFORMAT
#endif // DUMUX_IO_GRID_HH