gridmanager_yasp.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_IO_GRID_MANAGER_YASP_HH
#define DUMUX_IO_GRID_MANAGER_YASP_HH
 
#include <dune/common/math.hh>
 
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
 
#ifndef DUMUX_IO_GRID_MANAGER_BASE_HH
#include <dumux/io/grid/gridmanager_base.hh>
#endif
 
#include <dumux/common/gridcapabilities.hh>
 
namespace Dumux {
 
template<class Coordinates, int dim>
class GridManager<Dune::YaspGrid<dim, Coordinates>>
: public GridManagerBase<Dune::YaspGrid<dim, Coordinates>>
{
    using ct = typename Dune::YaspGrid<dim, Coordinates>::ctype;
    using GlobalPosition = Dune::FieldVector<ct, dim>;
public:
    using Grid = typename Dune::YaspGrid<dim, Coordinates>;
    using ParentType = GridManagerBase<Grid>;
 
    void init(const std::string& modelParamGroup = "")
    {
        // First try to create it from a DGF file in GridParameterGroup.File
        if (hasParamInGroup(modelParamGroup, "Grid.File"))
        {
            const auto filename = getParamFromGroup<std::string>(modelParamGroup, "Grid.File");
            const std::string extension = ParentType::getFileExtension(filename);
 
            if (extension != "dgf" && extension != "vti")
                DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " doesn't support grid files with extension: *."<< extension);
 
            // VTK file formats for image grids
            if (extension == "vti")
            {
#if DUMUX_HAVE_GRIDFORMAT
                VTKReader vtkReader(filename);
                VTKReader::Data cellData, pointData;
 
                const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
                std::shared_ptr<Detail::VTKReader::ImageGrid<ct, dim>> gridInfo
                    = vtkReader.readStructuredGrid<ct, dim>(cellData, pointData, verbose);
                const std::bitset<dim> periodic;
                const int overlap =  getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
                const auto upperRight = gridInfo->upperRight();
                const auto cells = gridInfo->cells();
                const auto lowerLeft = gridInfo->lowerLeft();
 
                if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>)
                    init_(upperRight, cells, modelParamGroup, overlap, periodic);
                else
                    init_(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic);
 
                ParentType::gridData_ = std::make_shared<GridData<Grid>>(
                    ParentType::gridPtr(), gridInfo, std::move(cellData), std::move(pointData)
                );
 
                ParentType::enableVtkData_ = true;
#else
                DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " doesn't support grid files with extension: *."<< extension);
#endif
            }
            else
            {
                ParentType::makeGridFromDgfFile(
                    getParamFromGroup<std::string>(modelParamGroup, "Grid.File")
                );
            }
 
            postProcessing_(modelParamGroup);
            return;
        }
 
        // Then look for the necessary keys to construct from the input file
        else if (hasParamInGroup(modelParamGroup, "Grid.UpperRight"))
        {
            // get the upper right corner coordinates
            const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight");
 
            // number of cells in each direction
            std::array<int, dim> cells; cells.fill(1);
            cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells);
 
            // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid)
            // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{});
            const std::bitset<dim> periodic;
 
            // get the overlap
            const int overlap =  getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
 
            if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>)
                init(upperRight, cells, modelParamGroup, overlap, periodic);
            else
            {
                const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0));
                init(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic);
            }
        }
 
        // Didn't find a way to construct the grid
        else
        {
            const auto prefix = modelParamGroup.empty() ? modelParamGroup : modelParamGroup + ".";
            DUNE_THROW(ParameterException, "Please supply one of the parameters "
                                           << prefix + "Grid.UpperRight"
                                           << ", or a grid file in " << prefix + "Grid.File");
 
        }
    }
 
    void init(const GlobalPosition& upperRight,
              const std::array<int, dim>& cells,
              const std::string& modelParamGroup = "",
              const int overlap = 1,
              const std::bitset<dim> periodic = std::bitset<dim>{})
    {
        init_(upperRight, cells, modelParamGroup, overlap, periodic);
        postProcessing_(modelParamGroup);
    }
 
    void init(const GlobalPosition& lowerLeft,
              const GlobalPosition& upperRight,
              const std::array<int, dim>& cells,
              const std::string& modelParamGroup = "",
              const int overlap = 1,
              const std::bitset<dim> periodic = std::bitset<dim>{})
    {
        init_(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic);
        postProcessing_(modelParamGroup);
    }
 
private:
    void init_(const GlobalPosition& upperRight,
               const std::array<int, dim>& cells,
               const std::string& modelParamGroup = "",
               const int overlap = 1,
               const std::bitset<dim> periodic = std::bitset<dim>{})
    {
        static_assert(std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>,
                      "Use init function taking lowerLeft as argument when working with EquidistantOffsetCoordinates");
 
        if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning"))
        {
            // construct using default load balancing
            ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap);
        }
        else
        {
            // construct using user defined partitioning
            const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
            Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning);
            ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb);
        }
    }
 
    void init_(const GlobalPosition& lowerLeft,
               const GlobalPosition& upperRight,
               const std::array<int, dim>& cells,
               const std::string& modelParamGroup = "",
               const int overlap = 1,
               const std::bitset<dim> periodic = std::bitset<dim>{})
    {
        static_assert(std::is_same_v<Dune::EquidistantOffsetCoordinates<ct, dim>, Coordinates>,
                      "LowerLeft can only be specified with EquidistantOffsetCoordinates");
 
        if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning"))
        {
            // construct using default load balancing
            ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap);
        }
        else
        {
            // construct using user defined partitioning
            const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
            Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning);
            ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb);
        }
    }
 
    void postProcessing_(const std::string& modelParamGroup)
    {
        // Check if should refine the grid
        const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true);
        ParentType::grid().refineOptions(keepPhysicalOverlap);
        ParentType::maybeRefineGrid(modelParamGroup);
        ParentType::loadBalance();
    }
};
 
template<class ctype, int dim>
class GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >>
: public GridManagerBase<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > >
{
public:
    using Grid = typename Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >;
    using ParentType = GridManagerBase<Grid>;
 
    void init(const std::string& modelParamGroup = "")
    {
        // Only construction from the input file is possible
        // Look for the necessary keys to construct from the input file
        // The positions
        std::array<std::vector<ctype>, dim> positions;
        for (int i = 0; i < dim; ++i)
            positions[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Positions" + std::to_string(i));
 
        // the number of cells (has a default)
        std::array<std::vector<int>, dim> cells;
        for (int i = 0; i < dim; ++i)
        {
            cells[i].resize(positions[i].size()-1, 1.0);
            cells[i] = getParamFromGroup<std::vector<int>>(modelParamGroup, "Grid.Cells" + std::to_string(i), cells[i]);
        }
 
        // grading factor (has a default)
        std::array<std::vector<ctype>, dim> grading;
        for (int i = 0; i < dim; ++i)
        {
            grading[i].resize(positions[i].size()-1, 1.0);
            grading[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]);
        }
 
        // call the generic function
        init(positions, cells, grading, modelParamGroup);
    }
 
    void init(const std::array<std::vector<ctype>, dim>& positions,
              const std::array<std::vector<int>, dim>& cells,
              const std::array<std::vector<ctype>, dim>& grading,
              const std::string& modelParamGroup = "")
    {
        // Additional parameters (they have a default)
        const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
        const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
        // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid)
        // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{});
        const std::bitset<dim> periodic;
 
        // Some sanity checks
        for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
        {
            if (cells[dimIdx].size() + 1 != positions[dimIdx].size())
            {
                DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Cells\" and \"Positions\" arrays");
            }
            if (grading[dimIdx].size() + 1 != positions[dimIdx].size())
            {
                DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Grading\" and \"Positions\" arrays");
            }
            ctype temp = std::numeric_limits<ctype>::lowest();
            for (unsigned int posIdx = 0; posIdx < positions[dimIdx].size(); ++posIdx)
            {
                if (temp > positions[dimIdx][posIdx])
                {
                    DUNE_THROW(Dune::RangeError, "Make sure to specify a monotone increasing \"Positions\" array");
                }
                temp = positions[dimIdx][posIdx];
            }
        }
 
        const auto globalPositions = computeGlobalPositions_(positions, cells, grading, verbose);
 
        // make the grid
        if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning"))
        {
            // construct using default load balancing
            ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap);
        }
        else
        {
            // construct using user defined partitioning
            const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
            Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning);
            ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap, typename Grid::Communication(), &lb);
        }
 
        postProcessing_(modelParamGroup);
    }
 
private:
    void postProcessing_(const std::string& modelParamGroup)
    {
        // Check if should refine the grid
        const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true);
        ParentType::grid().refineOptions(keepPhysicalOverlap);
        ParentType::maybeRefineGrid(modelParamGroup);
        ParentType::loadBalance();
    }
 
    std::array<std::vector<ctype>, dim>
    computeGlobalPositions_(const std::array<std::vector<ctype>, dim>& positions,
                            const std::array<std::vector<int>, dim>& cells,
                            const std::array<std::vector<ctype>, dim>& grading,
                            bool verbose = false)
    {
        std::array<std::vector<ctype>, dim> globalPositions;
        using std::pow;
        using Dune::power;
        for (int dimIdx = 0; dimIdx < dim; dimIdx++)
        {
            for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
            {
                ctype lower = positions[dimIdx][zoneIdx];
                ctype upper = positions[dimIdx][zoneIdx+1];
                int numCells = cells[dimIdx][zoneIdx];
                ctype gradingFactor = grading[dimIdx][zoneIdx];
                ctype length = upper - lower;
                ctype height = 1.0;
                bool increasingCellSize = false;
 
                if (verbose)
                {
                    std::cout << "dim " << dimIdx
                              << " lower "  << lower
                              << " upper "  << upper
                              << " numCells "  << numCells
                              << " grading "  << gradingFactor;
                }
 
                if (gradingFactor > 1.0)
                {
                    increasingCellSize = true;
                }
 
                // take absolute values and reverse cell size increment to achieve
                // reverse behavior for negative values
                if (gradingFactor < 0.0)
                {
                    using std::abs;
                    gradingFactor = abs(gradingFactor);
                    if (gradingFactor < 1.0)
                    {
                        increasingCellSize = true;
                    }
                }
 
                // if the grading factor is exactly 1.0 do equal spacing
                if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7)
                {
                    height = 1.0 / numCells;
                    if (verbose)
                    {
                        std::cout << " -> h "  << height * length << std::endl;
                    }
                }
                // if grading factor is not 1.0, do power law spacing
                else
                {
                    height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells));
 
                    if (verbose)
                    {
                        std::cout << " -> grading_eff "  << gradingFactor
                                  << " h_min "  << height * power(gradingFactor, 0) * length
                                  << " h_max "  << height * power(gradingFactor, numCells-1) * length
                                  << std::endl;
                    }
                }
 
                std::vector<ctype> localPositions;
                localPositions.push_back(0);
                for (int i = 0; i < numCells-1; i++)
                {
                    ctype hI = height;
                    if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7))
                    {
                        if (increasingCellSize)
                        {
                            hI *= power(gradingFactor, i);
                        }
                        else
                        {
                            hI *= power(gradingFactor, numCells-i-1);
                        }
                    }
                    localPositions.push_back(localPositions[i] + hI);
                }
 
                for (int i = 0; i < localPositions.size(); i++)
                {
                    localPositions[i] *= length;
                    localPositions[i] += lower;
                }
 
                for (unsigned int i = 0; i < localPositions.size(); ++i)
                {
                    globalPositions[dimIdx].push_back(localPositions[i]);
                }
            }
            globalPositions[dimIdx].push_back(positions[dimIdx].back());
        }
 
        return globalPositions;
    }
};
 
namespace Grid::Capabilities {
 
// To the best of our knowledge YaspGrid is view thread-safe
// This specialization can be removed after we depend on Dune release 2.9 and this is guaranteed by YaspGrid itself
template<class Coordinates, int dim>
struct MultithreadingSupported<Dune::YaspGrid<dim, Coordinates>>
{
    template<class GV>
    static bool eval(const GV& gv) // default is independent of the grid view
    { return true; }
};
 
} // end namespace Grid::Capabilities
 
} // end namespace Dumux
 
#endif