|
| 1 | +/* |
| 2 | + Copyright 2026 SINTEF Digital, Mathematics and Cybernetics. |
| 3 | +
|
| 4 | + This file is part of the Open Porous Media project (OPM). |
| 5 | +
|
| 6 | + OPM is free software: you can redistribute it and/or modify |
| 7 | + it under the terms of the GNU General Public License as published by |
| 8 | + the Free Software Foundation, either version 3 of the License, or |
| 9 | + (at your option) any later version. |
| 10 | +
|
| 11 | + OPM is distributed in the hope that it will be useful, |
| 12 | + but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 14 | + GNU General Public License for more details. |
| 15 | +
|
| 16 | + You should have received a copy of the GNU General Public License |
| 17 | + along with OPM. If not, see <http://www.gnu.org/licenses/>. |
| 18 | +*/ |
| 19 | + |
| 20 | +#include <config.h> |
| 21 | + |
| 22 | +#define NVERBOSE |
| 23 | + |
| 24 | +#define BOOST_TEST_MODULE PartitionOfUnityTest |
| 25 | +#define BOOST_TEST_NO_MAIN |
| 26 | +#include <boost/test/unit_test.hpp> |
| 27 | + |
| 28 | +#include <dune/common/parallel/mpihelper.hh> |
| 29 | + |
| 30 | +#include <opm/grid/CpGrid.hpp> |
| 31 | + |
| 32 | +#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp> |
| 33 | + |
| 34 | +#include <array> |
| 35 | +#include <vector> |
| 36 | + |
| 37 | +namespace |
| 38 | +{ |
| 39 | + |
| 40 | +/// Build a faulted 4x4x2 corner-point grid: the right half (i >= 2) is |
| 41 | +/// shifted down by half a cell. Face processing along the fault then |
| 42 | +/// produces nodes that are NOT among the eight canonical corners of the |
| 43 | +/// neighbouring cells (hanging nodes) -- the situation in which point |
| 44 | +/// (codim 3) communication interfaces built from cell_to_point_ miss |
| 45 | +/// nodes entirely. |
| 46 | +Opm::EclipseGrid faultedGrid(const int nx, const int ny, const int nz) |
| 47 | +{ |
| 48 | + const double shift = 0.5; |
| 49 | + |
| 50 | + // Straight vertical pillars on the unit lattice. |
| 51 | + std::vector<double> coord; |
| 52 | + coord.reserve((nx + 1) * (ny + 1) * 6); |
| 53 | + for (int j = 0; j <= ny; ++j) { |
| 54 | + for (int i = 0; i <= nx; ++i) { |
| 55 | + coord.push_back(i); // x top |
| 56 | + coord.push_back(j); // y top |
| 57 | + coord.push_back(0.0); // z top |
| 58 | + coord.push_back(i); // x bottom |
| 59 | + coord.push_back(j); // y bottom |
| 60 | + coord.push_back(nz + 1.0); // z bottom |
| 61 | + } |
| 62 | + } |
| 63 | + |
| 64 | + // ZCORN in Eclipse ordering: per layer, top surface then bottom surface, |
| 65 | + // each surface row-wise with two entries per cell in each direction. |
| 66 | + std::vector<double> zcorn; |
| 67 | + zcorn.reserve(8 * nx * ny * nz); |
| 68 | + for (int k = 0; k < nz; ++k) { |
| 69 | + for (int face = 0; face < 2; ++face) { // 0 = top, 1 = bottom |
| 70 | + for (int j = 0; j < ny; ++j) { |
| 71 | + for (int jj = 0; jj < 2; ++jj) { |
| 72 | + (void) jj; |
| 73 | + for (int i = 0; i < nx; ++i) { |
| 74 | + const double dz = (i >= nx / 2) ? shift : 0.0; |
| 75 | + const double z = k + face + dz; |
| 76 | + zcorn.push_back(z); |
| 77 | + zcorn.push_back(z); |
| 78 | + } |
| 79 | + } |
| 80 | + } |
| 81 | + } |
| 82 | + } |
| 83 | + |
| 84 | + return Opm::EclipseGrid(std::array<int, 3>{nx, ny, nz}, coord, zcorn); |
| 85 | +} |
| 86 | + |
| 87 | +/// Codim-3 data handle: senders write their rank for every vertex in the |
| 88 | +/// communication interface; receivers count the messages per vertex. |
| 89 | +class VertexReachedHandle |
| 90 | +{ |
| 91 | +public: |
| 92 | + using DataType = int; |
| 93 | + |
| 94 | + VertexReachedHandle(const Dune::CpGrid& grid, |
| 95 | + std::vector<int>& received, |
| 96 | + const int rank) |
| 97 | + : grid_(grid), received_(received), rank_(rank) |
| 98 | + {} |
| 99 | + |
| 100 | + bool fixedSize(int /*dim*/, int /*codim*/) |
| 101 | + { |
| 102 | + return true; |
| 103 | + } |
| 104 | + bool contains(int dim, int codim) |
| 105 | + { |
| 106 | + return dim == 3 && codim == 3; |
| 107 | + } |
| 108 | + template<class T> |
| 109 | + std::size_t size(const T&) |
| 110 | + { |
| 111 | + return 1; |
| 112 | + } |
| 113 | + template<class B, class T> |
| 114 | + void gather(B& buffer, const T&) |
| 115 | + { |
| 116 | + buffer.write(rank_); |
| 117 | + } |
| 118 | + template<class B, class T> |
| 119 | + void scatter(B& buffer, const T& t, std::size_t) |
| 120 | + { |
| 121 | + int sender = -1; |
| 122 | + buffer.read(sender); |
| 123 | + ++received_[grid_.leafIndexSet().index(t)]; |
| 124 | + } |
| 125 | + |
| 126 | +private: |
| 127 | + const Dune::CpGrid& grid_; |
| 128 | + std::vector<int>& received_; |
| 129 | + int rank_; |
| 130 | +}; |
| 131 | + |
| 132 | +} // anonymous namespace |
| 133 | + |
| 134 | +// Every vertex on the processor boundary (partition type != interior) must be |
| 135 | +// reached by codim-3 communication. On a faulted grid the fault-face hanging |
| 136 | +// nodes are not canonical corners of the neighbouring cells; interfaces built |
| 137 | +// from the canonical corners only never include them, so they receive no |
| 138 | +// message and their ownership cannot be made consistent across ranks (the |
| 139 | +// "Owner is not a partition of unity" failure seen in vertex-based mechanics). |
| 140 | +BOOST_AUTO_TEST_CASE(FaultedGridVerticesReached) |
| 141 | +{ |
| 142 | + const auto& helper = Dune::MPIHelper::instance( |
| 143 | + boost::unit_test::framework::master_test_suite().argc, |
| 144 | + boost::unit_test::framework::master_test_suite().argv); |
| 145 | + const int size = helper.size(); |
| 146 | + const int rank = helper.rank(); |
| 147 | + |
| 148 | + if (size < 2) { |
| 149 | + return; // needs at least two processes to have processor boundaries |
| 150 | + } |
| 151 | + |
| 152 | + const int nx = 4, ny = 4, nz = 2; |
| 153 | + const auto eclGrid = faultedGrid(nx, ny, nz); |
| 154 | + |
| 155 | + Dune::CpGrid grid; |
| 156 | + grid.processEclipseFormat(&eclGrid, nullptr, |
| 157 | + /* periodic_extension = */ false, |
| 158 | + /* turn_normals = */ false, |
| 159 | + /* clip_z = */ false, |
| 160 | + /* pinchActive = */ false, |
| 161 | + /* edge_conformal = */ false); |
| 162 | + |
| 163 | + // Partition in vertical slabs of columns so the fault plane (i == nx/2) |
| 164 | + // coincides with a processor boundary. |
| 165 | + std::vector<int> parts(grid.size(0)); |
| 166 | + const auto& gv = grid.leafGridView(); |
| 167 | + for (const auto& element : elements(gv)) { |
| 168 | + const int cell = gv.indexSet().index(element); |
| 169 | + const auto center = element.geometry().center(); |
| 170 | + const int i = static_cast<int>(center[0]); // unit lattice |
| 171 | + const int j = static_cast<int>(center[1]); |
| 172 | + if (size >= 4) { |
| 173 | + parts[cell] = (i >= nx / 2) + 2 * (j >= ny / 2); |
| 174 | + } else { |
| 175 | + parts[cell] = (i >= nx / 2); |
| 176 | + } |
| 177 | + } |
| 178 | + |
| 179 | + grid.loadBalance(parts, /* ownersFirst = */ false, |
| 180 | + /* addCornerCells = */ true, /* overlapLayers = */ 1); |
| 181 | + |
| 182 | + const int numVertices = grid.leafIndexSet().size(3); |
| 183 | + std::vector<int> received(numVertices, 0); |
| 184 | + |
| 185 | + VertexReachedHandle handle(grid, received, rank); |
| 186 | + grid.communicate(handle, Dune::All_All_Interface, |
| 187 | + Dune::ForwardCommunication); |
| 188 | + |
| 189 | + // Every non-interior vertex sits on a processor boundary and is shared |
| 190 | + // with at least one other rank, so it must have received at least one |
| 191 | + // message. Hanging nodes on the fault violate this if the interface was |
| 192 | + // built from the canonical corners only. |
| 193 | + int unreached = 0; |
| 194 | + for (const auto& vertex : vertices(grid.leafGridView())) { |
| 195 | + const int idx = grid.leafIndexSet().index(vertex); |
| 196 | + if (vertex.partitionType() != Dune::InteriorEntity |
| 197 | + && received[idx] == 0) { |
| 198 | + ++unreached; |
| 199 | + } |
| 200 | + } |
| 201 | + |
| 202 | + BOOST_CHECK_MESSAGE(unreached == 0, |
| 203 | + "rank " << rank << ": " << unreached |
| 204 | + << " non-interior vertices were not reached by codim-3" |
| 205 | + " communication (hanging nodes missing from the point" |
| 206 | + " interface)"); |
| 207 | +} |
| 208 | + |
| 209 | +bool |
| 210 | +init_unit_test_func() |
| 211 | +{ |
| 212 | + return true; |
| 213 | +} |
| 214 | + |
| 215 | +int main(int argc, char** argv) |
| 216 | +{ |
| 217 | + Dune::MPIHelper::instance(argc, argv); |
| 218 | + |
| 219 | + return boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv); |
| 220 | +} |
0 commit comments