| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup InputOutput | ||
| 10 | * \brief Grid manager specialization for YaspGrid | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_IO_GRID_MANAGER_YASP_HH | ||
| 13 | #define DUMUX_IO_GRID_MANAGER_YASP_HH | ||
| 14 | |||
| 15 | #include <dune/common/math.hh> | ||
| 16 | |||
| 17 | #include <dune/grid/yaspgrid.hh> | ||
| 18 | #include <dune/grid/io/file/dgfparser/dgfyasp.hh> | ||
| 19 | |||
| 20 | #ifndef DUMUX_IO_GRID_MANAGER_BASE_HH | ||
| 21 | #include <dumux/io/grid/gridmanager_base.hh> | ||
| 22 | #endif | ||
| 23 | |||
| 24 | #include <dumux/common/gridcapabilities.hh> | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup InputOutput | ||
| 30 | * \brief Provides a grid manager for YaspGrids | ||
| 31 | * from information in the input file | ||
| 32 | * | ||
| 33 | * All keys are expected to be in group GridParameterGroup. | ||
| 34 | * The following keys are recognized: | ||
| 35 | * - File : a DGF file to load the coarse grid from | ||
| 36 | * - LowerLeft : lower left corner of the domain (only for EquidistantOffsetCoordinates, otherwise 0-origin) | ||
| 37 | * - UpperRight : upper right corner of the domain | ||
| 38 | * - Cells : the number of cells in each direction | ||
| 39 | * - Periodic : true or false for each direction | ||
| 40 | * - Overlap : overlap size in cells | ||
| 41 | * - Partitioning : a non-standard load-balancing, number of processors per direction | ||
| 42 | * - KeepPyhsicalOverlap : whether to keep the physical overlap | ||
| 43 | * in physical size or in number of cells upon refinement | ||
| 44 | * - Refinement : the number of global refines to apply initially. | ||
| 45 | * | ||
| 46 | */ | ||
| 47 | template<class Coordinates, int dim> | ||
| 48 |
3/5✓ Branch 1 taken 29 times.
✓ Branch 2 taken 31 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
571 | class GridManager<Dune::YaspGrid<dim, Coordinates>> |
| 49 | : public GridManagerBase<Dune::YaspGrid<dim, Coordinates>> | ||
| 50 | { | ||
| 51 | using ct = typename Dune::YaspGrid<dim, Coordinates>::ctype; | ||
| 52 | using GlobalPosition = Dune::FieldVector<ct, dim>; | ||
| 53 | public: | ||
| 54 | using Grid = typename Dune::YaspGrid<dim, Coordinates>; | ||
| 55 | using ParentType = GridManagerBase<Grid>; | ||
| 56 | |||
| 57 | /*! | ||
| 58 | * \brief Make the grid. This is implemented by specializations of this method. | ||
| 59 | */ | ||
| 60 | 375 | void init(const std::string& modelParamGroup = "") | |
| 61 | { | ||
| 62 | // First try to create it from a DGF file in GridParameterGroup.File | ||
| 63 |
6/14✓ Branch 1 taken 313 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 313 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 313 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 313 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 313 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 313 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
1125 | if (hasParamInGroup(modelParamGroup, "Grid.File")) |
| 64 | { | ||
| 65 | ✗ | ParentType::makeGridFromDgfFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File")); | |
| 66 | ✗ | postProcessing_(modelParamGroup); | |
| 67 | ✗ | return; | |
| 68 | } | ||
| 69 | |||
| 70 | // Then look for the necessary keys to construct from the input file | ||
| 71 |
6/14✓ Branch 1 taken 313 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 313 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 313 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 313 times.
✓ Branch 11 taken 313 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 313 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
1125 | else if (hasParamInGroup(modelParamGroup, "Grid.UpperRight")) |
| 72 | { | ||
| 73 | // get the upper right corner coordinates | ||
| 74 | 375 | const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight"); | |
| 75 | |||
| 76 | // number of cells in each direction | ||
| 77 | 375 | std::array<int, dim> cells; cells.fill(1); | |
| 78 | 375 | cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells); | |
| 79 | |||
| 80 | // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid) | ||
| 81 | // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{}); | ||
| 82 | 375 | const std::bitset<dim> periodic; | |
| 83 | |||
| 84 | // get the overlap | ||
| 85 | 375 | const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); | |
| 86 | |||
| 87 | if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>) | ||
| 88 | 253 | init(upperRight, cells, modelParamGroup, overlap, periodic); | |
| 89 | else | ||
| 90 | { | ||
| 91 | 244 | const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0)); | |
| 92 | 122 | init(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic); | |
| 93 | } | ||
| 94 | } | ||
| 95 | |||
| 96 | // Didn't find a way to construct the grid | ||
| 97 | else | ||
| 98 | { | ||
| 99 | ✗ | const auto prefix = modelParamGroup.empty() ? modelParamGroup : modelParamGroup + "."; | |
| 100 | ✗ | DUNE_THROW(ParameterException, "Please supply one of the parameters " | |
| 101 | << prefix + "Grid.UpperRight" | ||
| 102 | << ", or a grid file in " << prefix + "Grid.File"); | ||
| 103 | |||
| 104 | } | ||
| 105 | } | ||
| 106 | |||
| 107 | /*! | ||
| 108 | * \brief Make the grid using input data not read from the input file. | ||
| 109 | * \note Use this function for EquidistantCoordinates where lowerLeft is always zero. | ||
| 110 | */ | ||
| 111 | 224 | void init(const GlobalPosition& upperRight, | |
| 112 | const std::array<int, dim>& cells, | ||
| 113 | const std::string& modelParamGroup = "", | ||
| 114 | const int overlap = 1, | ||
| 115 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 116 | { | ||
| 117 | static_assert(std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>, | ||
| 118 | "Use init function taking lowerLeft as argument when working with EquidistantOffsetCoordinates"); | ||
| 119 | |||
| 120 |
6/14✓ Branch 1 taken 220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 220 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 220 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 220 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 220 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 220 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
1120 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 121 | { | ||
| 122 | // construct using default load balancing | ||
| 123 |
3/8✓ Branch 2 taken 220 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 220 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 220 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
224 | ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap); |
| 124 | } | ||
| 125 | else | ||
| 126 | { | ||
| 127 | // construct using user defined partitioning | ||
| 128 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 129 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 130 | ✗ | ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb); | |
| 131 | } | ||
| 132 | |||
| 133 | 224 | postProcessing_(modelParamGroup); | |
| 134 | 224 | } | |
| 135 | |||
| 136 | /*! | ||
| 137 | * \brief Make the grid using input data not read from the input file. | ||
| 138 | * \note Use this function for EquidistantOffsetCoordinates. | ||
| 139 | */ | ||
| 140 | 95 | void init(const GlobalPosition& lowerLeft, | |
| 141 | const GlobalPosition& upperRight, | ||
| 142 | const std::array<int, dim>& cells, | ||
| 143 | const std::string& modelParamGroup = "", | ||
| 144 | const int overlap = 1, | ||
| 145 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 146 | { | ||
| 147 | static_assert(std::is_same_v<Dune::EquidistantOffsetCoordinates<ct, dim>, Coordinates>, | ||
| 148 | "LowerLeft can only be specified with EquidistantOffsetCoordinates"); | ||
| 149 | |||
| 150 |
6/14✓ Branch 1 taken 95 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 95 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 95 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 95 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 95 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 95 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
475 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 151 | { | ||
| 152 | // construct using default load balancing | ||
| 153 |
3/8✓ Branch 2 taken 95 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 95 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 95 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
95 | ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap); |
| 154 | } | ||
| 155 | else | ||
| 156 | { | ||
| 157 | // construct using user defined partitioning | ||
| 158 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 159 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 160 | ✗ | ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb); | |
| 161 | } | ||
| 162 | |||
| 163 | 95 | postProcessing_(modelParamGroup); | |
| 164 | 95 | } | |
| 165 | |||
| 166 | private: | ||
| 167 | |||
| 168 | /*! | ||
| 169 | * \brief Postprocessing for YaspGrid | ||
| 170 | */ | ||
| 171 | 377 | void postProcessing_(const std::string& modelParamGroup) | |
| 172 | { | ||
| 173 | // Check if should refine the grid | ||
| 174 | 377 | const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true); | |
| 175 | 377 | ParentType::grid().refineOptions(keepPhysicalOverlap); | |
| 176 | 377 | ParentType::maybeRefineGrid(modelParamGroup); | |
| 177 | 377 | ParentType::loadBalance(); | |
| 178 | 377 | } | |
| 179 | }; | ||
| 180 | |||
| 181 | /*! | ||
| 182 | * \ingroup InputOutput | ||
| 183 | * \brief Provides a grid manager for YaspGrids with different zones and grading | ||
| 184 | * | ||
| 185 | * All keys are expected to be in group GridParameterGroup. | ||
| 186 | * The following keys are recognized: | ||
| 187 | * - Positions0 : position array for x-coordinate | ||
| 188 | * - Positions1 : position array for y-coordinate | ||
| 189 | * - Positions2 : position array for z-coordinate | ||
| 190 | * - Cells0 : number of cells array for x-coordinate | ||
| 191 | * - Cells1 : number of cells array for y-coordinate | ||
| 192 | * - Cells2 : number of cells array for z-coordinate | ||
| 193 | * - Grading0 : grading factor array for x-coordinate | ||
| 194 | * - Grading1 : grading factor array for y-coordinate | ||
| 195 | * - Grading2 : grading factor array for z-coordinate | ||
| 196 | * - Verbosity : whether the grid construction should output to standard out | ||
| 197 | * - Periodic : true or false for each direction | ||
| 198 | * - Overlap : overlap size in cells | ||
| 199 | * - Partitioning : a non-standard load-balancing, number of processors per direction | ||
| 200 | * - KeepPyhsicalOverlap : whether to keep the physical overlap | ||
| 201 | * in physical size or in number of cells upon refinement | ||
| 202 | * - Refinement : the number of global refines to apply initially. | ||
| 203 | * | ||
| 204 | * The grading factor \f$ g \f$ specifies the ratio between the next and the current cell size: | ||
| 205 | * \f$ g = \frac{h_{i+1}}{h_i} \f$. | ||
| 206 | * Negative grading factors are converted to | ||
| 207 | * \f$ g = -\frac{1}{g_\textrm{negative}} \f$ | ||
| 208 | * to avoid issues with imprecise fraction numbers. | ||
| 209 | */ | ||
| 210 | template<class ctype, int dim> | ||
| 211 |
4/7✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 8 times.
✗ Branch 14 not taken.
|
165 | class GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >> |
| 212 | : public GridManagerBase<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > > | ||
| 213 | { | ||
| 214 | public: | ||
| 215 | using Grid = typename Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >; | ||
| 216 | using ParentType = GridManagerBase<Grid>; | ||
| 217 | |||
| 218 | /*! | ||
| 219 | * \brief Make the grid. This is implemented by specializations of this method. | ||
| 220 | */ | ||
| 221 | 66 | void init(const std::string& modelParamGroup = "") | |
| 222 | { | ||
| 223 | // Only construction from the input file is possible | ||
| 224 | // Look for the necessary keys to construct from the input file | ||
| 225 | // The positions | ||
| 226 | 66 | std::array<std::vector<ctype>, dim> positions; | |
| 227 |
2/2✓ Branch 0 taken 66 times.
✓ Branch 1 taken 121 times.
|
187 | for (int i = 0; i < dim; ++i) |
| 228 |
6/18✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 121 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 121 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 121 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
363 | positions[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Positions" + std::to_string(i)); |
| 229 | |||
| 230 | // the number of cells (has a default) | ||
| 231 | 66 | std::array<std::vector<int>, dim> cells; | |
| 232 |
2/2✓ Branch 0 taken 66 times.
✓ Branch 1 taken 121 times.
|
187 | for (int i = 0; i < dim; ++i) |
| 233 | { | ||
| 234 |
4/8✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 121 times.
✗ Branch 11 not taken.
|
484 | cells[i].resize(positions[i].size()-1, 1.0); |
| 235 |
7/18✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 121 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 121 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 121 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 121 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
484 | cells[i] = getParamFromGroup<std::vector<int>>(modelParamGroup, "Grid.Cells" + std::to_string(i), cells[i]); |
| 236 | } | ||
| 237 | |||
| 238 | // grading factor (has a default) | ||
| 239 | 66 | std::array<std::vector<ctype>, dim> grading; | |
| 240 |
2/2✓ Branch 0 taken 121 times.
✓ Branch 1 taken 66 times.
|
187 | for (int i = 0; i < dim; ++i) |
| 241 | { | ||
| 242 |
4/8✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 121 times.
✗ Branch 11 not taken.
|
484 | grading[i].resize(positions[i].size()-1, 1.0); |
| 243 |
7/18✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 121 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 121 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 121 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 121 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
484 | grading[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]); |
| 244 | } | ||
| 245 | |||
| 246 | // call the generic function | ||
| 247 |
1/2✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
|
66 | init(positions, cells, grading, modelParamGroup); |
| 248 | 66 | } | |
| 249 | |||
| 250 | /*! | ||
| 251 | * \brief Make the grid using input data not read from the input file. | ||
| 252 | */ | ||
| 253 | 96 | void init(const std::array<std::vector<ctype>, dim>& positions, | |
| 254 | const std::array<std::vector<int>, dim>& cells, | ||
| 255 | const std::array<std::vector<ctype>, dim>& grading, | ||
| 256 | const std::string& modelParamGroup = "") | ||
| 257 | { | ||
| 258 | // Additional parameters (they have a default) | ||
| 259 | 96 | const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); | |
| 260 | 96 | const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false); | |
| 261 | // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid) | ||
| 262 | // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{}); | ||
| 263 | 96 | const std::bitset<dim> periodic; | |
| 264 | |||
| 265 | // Some sanity checks | ||
| 266 |
2/2✓ Branch 0 taken 170 times.
✓ Branch 1 taken 86 times.
|
288 | for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) |
| 267 | { | ||
| 268 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 170 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 170 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 170 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 170 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 170 times.
|
960 | if (cells[dimIdx].size() + 1 != positions[dimIdx].size()) |
| 269 | { | ||
| 270 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Cells\" and \"Positions\" arrays"); | |
| 271 | } | ||
| 272 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 170 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 170 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 170 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 170 times.
|
768 | if (grading[dimIdx].size() + 1 != positions[dimIdx].size()) |
| 273 | { | ||
| 274 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Grading\" and \"Positions\" arrays"); | |
| 275 | } | ||
| 276 | ctype temp = std::numeric_limits<ctype>::lowest(); | ||
| 277 |
4/4✓ Branch 0 taken 355 times.
✓ Branch 1 taken 181 times.
✓ Branch 2 taken 355 times.
✓ Branch 3 taken 181 times.
|
1204 | for (unsigned int posIdx = 0; posIdx < positions[dimIdx].size(); ++posIdx) |
| 278 | { | ||
| 279 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 366 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 366 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 366 times.
|
1230 | if (temp > positions[dimIdx][posIdx]) |
| 280 | { | ||
| 281 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify a monotone increasing \"Positions\" array"); | |
| 282 | } | ||
| 283 | 820 | temp = positions[dimIdx][posIdx]; | |
| 284 | } | ||
| 285 | } | ||
| 286 | |||
| 287 | 96 | const auto globalPositions = computeGlobalPositions_(positions, cells, grading, verbose); | |
| 288 | |||
| 289 | // make the grid | ||
| 290 |
6/14✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 86 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 86 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
480 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 291 | { | ||
| 292 | // construct using default load balancing | ||
| 293 |
4/12✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
96 | ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap); |
| 294 | } | ||
| 295 | else | ||
| 296 | { | ||
| 297 | // construct using user defined partitioning | ||
| 298 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 299 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 300 | ✗ | ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap, typename Grid::Communication(), &lb); | |
| 301 | } | ||
| 302 | |||
| 303 |
1/2✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
|
96 | postProcessing_(modelParamGroup); |
| 304 | 96 | } | |
| 305 | |||
| 306 | private: | ||
| 307 | /*! | ||
| 308 | * \brief Postprocessing for YaspGrid | ||
| 309 | */ | ||
| 310 | 96 | void postProcessing_(const std::string& modelParamGroup) | |
| 311 | { | ||
| 312 | // Check if should refine the grid | ||
| 313 | 96 | const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true); | |
| 314 | 96 | ParentType::grid().refineOptions(keepPhysicalOverlap); | |
| 315 | 96 | ParentType::maybeRefineGrid(modelParamGroup); | |
| 316 | 96 | ParentType::loadBalance(); | |
| 317 | 96 | } | |
| 318 | |||
| 319 | //! Compute the global position tensor grid from the given positions, cells, and grading factors | ||
| 320 | std::array<std::vector<ctype>, dim> | ||
| 321 | 86 | computeGlobalPositions_(const std::array<std::vector<ctype>, dim>& positions, | |
| 322 | const std::array<std::vector<int>, dim>& cells, | ||
| 323 | const std::array<std::vector<ctype>, dim>& grading, | ||
| 324 | bool verbose = false) | ||
| 325 | { | ||
| 326 | 86 | std::array<std::vector<ctype>, dim> globalPositions; | |
| 327 | using std::pow; | ||
| 328 | using Dune::power; | ||
| 329 |
2/2✓ Branch 0 taken 148 times.
✓ Branch 1 taken 76 times.
|
256 | for (int dimIdx = 0; dimIdx < dim; dimIdx++) |
| 330 | { | ||
| 331 |
6/6✓ Branch 0 taken 174 times.
✓ Branch 1 taken 148 times.
✓ Branch 2 taken 174 times.
✓ Branch 3 taken 148 times.
✓ Branch 4 taken 174 times.
✓ Branch 5 taken 148 times.
|
1098 | for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) |
| 332 | { | ||
| 333 |
4/4✓ Branch 0 taken 73 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 73 times.
✓ Branch 3 taken 101 times.
|
392 | ctype lower = positions[dimIdx][zoneIdx]; |
| 334 |
4/4✓ Branch 0 taken 73 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 73 times.
✓ Branch 3 taken 101 times.
|
392 | ctype upper = positions[dimIdx][zoneIdx+1]; |
| 335 |
4/4✓ Branch 0 taken 73 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 73 times.
✓ Branch 3 taken 101 times.
|
392 | int numCells = cells[dimIdx][zoneIdx]; |
| 336 |
4/4✓ Branch 0 taken 73 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 73 times.
✓ Branch 3 taken 101 times.
|
392 | ctype gradingFactor = grading[dimIdx][zoneIdx]; |
| 337 | 196 | ctype length = upper - lower; | |
| 338 | 196 | ctype height = 1.0; | |
| 339 | 196 | bool increasingCellSize = false; | |
| 340 | |||
| 341 |
2/2✓ Branch 0 taken 73 times.
✓ Branch 1 taken 101 times.
|
196 | if (verbose) |
| 342 | { | ||
| 343 | 73 | std::cout << "dim " << dimIdx | |
| 344 |
2/4✓ Branch 1 taken 73 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 73 times.
✗ Branch 6 not taken.
|
146 | << " lower " << lower |
| 345 |
1/2✓ Branch 2 taken 73 times.
✗ Branch 3 not taken.
|
146 | << " upper " << upper |
| 346 | 73 | << " numCells " << numCells | |
| 347 |
2/4✓ Branch 1 taken 73 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 73 times.
✗ Branch 6 not taken.
|
146 | << " grading " << gradingFactor; |
| 348 | } | ||
| 349 | |||
| 350 |
2/2✓ Branch 0 taken 27 times.
✓ Branch 1 taken 147 times.
|
196 | if (gradingFactor > 1.0) |
| 351 | { | ||
| 352 | 27 | increasingCellSize = true; | |
| 353 | } | ||
| 354 | |||
| 355 | // take absolute values and reverse cell size increment to achieve | ||
| 356 | // reverse behavior for negative values | ||
| 357 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 152 times.
|
196 | if (gradingFactor < 0.0) |
| 358 | { | ||
| 359 | using std::abs; | ||
| 360 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 17 times.
|
22 | gradingFactor = abs(gradingFactor); |
| 361 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 17 times.
|
22 | if (gradingFactor < 1.0) |
| 362 | { | ||
| 363 | 5 | increasingCellSize = true; | |
| 364 | } | ||
| 365 | } | ||
| 366 | |||
| 367 | // if the grading factor is exactly 1.0 do equal spacing | ||
| 368 |
4/4✓ Branch 0 taken 164 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 120 times.
✓ Branch 3 taken 44 times.
|
196 | if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7) |
| 369 | { | ||
| 370 | 142 | height = 1.0 / numCells; | |
| 371 |
2/2✓ Branch 0 taken 35 times.
✓ Branch 1 taken 85 times.
|
142 | if (verbose) |
| 372 | { | ||
| 373 |
2/4✓ Branch 2 taken 35 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 35 times.
✗ Branch 6 not taken.
|
70 | std::cout << " -> h " << height * length << std::endl; |
| 374 | } | ||
| 375 | } | ||
| 376 | // if grading factor is not 1.0, do power law spacing | ||
| 377 | else | ||
| 378 | { | ||
| 379 | 54 | height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells)); | |
| 380 | |||
| 381 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 16 times.
|
54 | if (verbose) |
| 382 | { | ||
| 383 |
1/2✓ Branch 2 taken 38 times.
✗ Branch 3 not taken.
|
76 | std::cout << " -> grading_eff " << gradingFactor |
| 384 |
2/4✓ Branch 2 taken 38 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 38 times.
✗ Branch 6 not taken.
|
114 | << " h_min " << height * power(gradingFactor, 0) * length |
| 385 |
1/2✓ Branch 2 taken 38 times.
✗ Branch 3 not taken.
|
114 | << " h_max " << height * power(gradingFactor, numCells-1) * length |
| 386 |
1/2✓ Branch 1 taken 38 times.
✗ Branch 2 not taken.
|
38 | << std::endl; |
| 387 | } | ||
| 388 | } | ||
| 389 | |||
| 390 |
2/4✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 174 times.
✗ Branch 4 not taken.
|
588 | std::vector<ctype> localPositions; |
| 391 |
1/2✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
|
196 | localPositions.push_back(0); |
| 392 |
2/2✓ Branch 0 taken 9861 times.
✓ Branch 1 taken 174 times.
|
10123 | for (int i = 0; i < numCells-1; i++) |
| 393 | { | ||
| 394 | 9927 | ctype hI = height; | |
| 395 |
4/4✓ Branch 0 taken 9331 times.
✓ Branch 1 taken 530 times.
✓ Branch 2 taken 5220 times.
✓ Branch 3 taken 4111 times.
|
9927 | if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7)) |
| 396 | { | ||
| 397 |
2/2✓ Branch 0 taken 4395 times.
✓ Branch 1 taken 1355 times.
|
5750 | if (increasingCellSize) |
| 398 | { | ||
| 399 | 8790 | hI *= power(gradingFactor, i); | |
| 400 | } | ||
| 401 | else | ||
| 402 | { | ||
| 403 | 2710 | hI *= power(gradingFactor, numCells-i-1); | |
| 404 | } | ||
| 405 | } | ||
| 406 |
2/6✓ Branch 1 taken 9861 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9861 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
19854 | localPositions.push_back(localPositions[i] + hI); |
| 407 | } | ||
| 408 | |||
| 409 |
4/4✓ Branch 0 taken 10035 times.
✓ Branch 1 taken 174 times.
✓ Branch 2 taken 10035 times.
✓ Branch 3 taken 174 times.
|
20638 | for (int i = 0; i < localPositions.size(); i++) |
| 410 | { | ||
| 411 | 10123 | localPositions[i] *= length; | |
| 412 | 10123 | localPositions[i] += lower; | |
| 413 | } | ||
| 414 | |||
| 415 |
4/4✓ Branch 0 taken 10035 times.
✓ Branch 1 taken 174 times.
✓ Branch 2 taken 10035 times.
✓ Branch 3 taken 174 times.
|
10515 | for (unsigned int i = 0; i < localPositions.size(); ++i) |
| 416 | { | ||
| 417 |
3/6✓ Branch 1 taken 10035 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10035 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10035 times.
✗ Branch 8 not taken.
|
30369 | globalPositions[dimIdx].push_back(localPositions[i]); |
| 418 | } | ||
| 419 | } | ||
| 420 |
4/8✓ Branch 1 taken 148 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 148 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 148 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 148 times.
✗ Branch 11 not taken.
|
680 | globalPositions[dimIdx].push_back(positions[dimIdx].back()); |
| 421 | } | ||
| 422 | |||
| 423 | 86 | return globalPositions; | |
| 424 | } | ||
| 425 | }; | ||
| 426 | |||
| 427 | namespace Grid::Capabilities { | ||
| 428 | |||
| 429 | // To the best of our knowledge YaspGrid is view thread-safe | ||
| 430 | // This specialization can be removed after we depend on Dune release 2.9 and this is guaranteed by YaspGrid itself | ||
| 431 | template<class Coordinates, int dim> | ||
| 432 | struct MultithreadingSupported<Dune::YaspGrid<dim, Coordinates>> | ||
| 433 | { | ||
| 434 | template<class GV> | ||
| 435 | ✗ | static bool eval(const GV& gv) // default is independent of the grid view | |
| 436 | ✗ | { return true; } | |
| 437 | }; | ||
| 438 | |||
| 439 | } // end namespace Grid::Capabilities | ||
| 440 | |||
| 441 | } // end namespace Dumux | ||
| 442 | |||
| 443 | #endif | ||
| 444 |