diff --git a/src/opm/output/eclipse/ActiveIndexByColumns.cpp b/src/opm/output/eclipse/ActiveIndexByColumns.cpp index b4354bdc7..fc554831b 100644 --- a/src/opm/output/eclipse/ActiveIndexByColumns.cpp +++ b/src/opm/output/eclipse/ActiveIndexByColumns.cpp @@ -26,11 +26,14 @@ #include #include #include +#include #include namespace { - std::size_t columnarGlobalIdx(const std::array& dims, - const std::array& ijk) + std::size_t columnarGlobalIdx(const std::array& dims, + const std::array& ijk, + const std::array::size_type outer, + const std::array::size_type middle) { // Linear index assuming C-like loop order // @@ -38,13 +41,38 @@ namespace { // for (j = 0 .. Ny - 1) // for (k = 0 .. Nz - 1) // - // as opposed to the usual Fortran-like loop order ("natural ordering") + // or, if Ny > Nx, the alternate loop order + // + // for (j = 0 .. Ny - 1) + // for (i = 0 .. Nx - 1) + // for (k = 0 .. Nz - 1) + // + // swapping the 'i' and 'j' loops. This is opposed to the usual, + // Fortran-like, loop order ("natural ordering") // // for (k = 0 .. Nz - 1) // for (j = 0 .. Ny - 1) // for (i = 0 .. Nx - 1) // - return ijk[2] + dims[2]*(ijk[1] + dims[1]*ijk[0]); + // that is used elsewhere. In other words, the inner loop is always + // across the model layers while we ensure that the 'outer' loop + // always iterates over an index range that is at least as large as + // that of the 'middle' loop. + + return ijk[2] + dims[2]*(ijk[middle] + dims[middle]*ijk[outer]); + } + + std::pair::size_type, std::array::size_type> + inferOuterLoopOrdering(const std::array& cartDims) + { + auto outer = std::array::size_type{0}; + auto middle = std::array::size_type{1}; + + if (cartDims[middle] > cartDims[outer]) { + std::swap(outer, middle); + } + + return { outer, middle }; } std::vector @@ -54,11 +82,13 @@ namespace { { auto colGlobIx = activeCells; + const auto [outer, middle] = inferOuterLoopOrdering(cartDims); + std::transform(colGlobIx.begin(), colGlobIx.end(), colGlobIx.begin(), - [&cartDims, &getIJK] + [&cartDims, &getIJK, outer, middle] (const std::size_t cell) { - return columnarGlobalIdx(cartDims, getIJK(cell)); + return columnarGlobalIdx(cartDims, getIJK(cell), outer, middle); }); return colGlobIx; diff --git a/tests/test_ActiveIndexByColumns.cpp b/tests/test_ActiveIndexByColumns.cpp index 8674a84a2..d5af75045 100644 --- a/tests/test_ActiveIndexByColumns.cpp +++ b/tests/test_ActiveIndexByColumns.cpp @@ -477,3 +477,414 @@ BOOST_AUTO_TEST_CASE(Cube_3x3x3_exclude_diagonals) } BOOST_AUTO_TEST_SUITE_END() // Grid_Based + +// ===================================================================== + +BOOST_AUTO_TEST_SUITE(Grid_Based_NY_Larger_Than_NX) + +namespace { + std::vector coord_2x3x4() + { + return { + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 2.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 2.0, 1.0, 0.0, 2.0, 1.0, 0.0, + 0.0, 2.0, 0.0, 0.0, 2.0, 0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.0, 2.0, 2.0, 0.0, 2.0, 2.0, 0.0, + 0.0, 3.0, 0.0, 0.0, 3.0, 0.0, 1.0, 3.0, 0.0, 1.0, 3.0, 0.0, 2.0, 3.0, 0.0, 2.0, 3.0, 0.0, + }; + } + + std::vector zcorn_2x3x4() + { + return { + // Top, layer 1 + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // 0.. 1 + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // 2.. 3 + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // 4.. 5 + + // Bottom, layer 1 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 0.. 1 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 2.. 3 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 4.. 5 + + // Top, layer 2 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 6.. 7 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 8.. 9 + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 10..11 + + // Bottom, layer 2 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 6.. 7 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 8.. 9 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 10..11 + + // Top, layer 3 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 12..13 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 14..15 + 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, // 16..17 + + // Bottom, layer 3 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 12..13 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 14..15 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 16..17 + + // Top, layer 4 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 18..19 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 20..21 + 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, // 22..23 + + // Bottom, layer 4 + 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, // 18..19 + 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, // 20..21 + 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, // 22..23 + }; + } + + std::vector actnum_2x3x4_exclude_alternate_centre() + { + return { + 1, 1, + 0, 1, + 1, 1, + + 1, 1, + 1, 0, + 1, 1, + + 1, 1, + 0, 1, + 1, 1, + + 1, 1, + 1, 0, + 1, 1, + }; + } + + std::vector actnum_2x3x4_exclude_centre_column() + { + return { + 1, 1, + 1, 0, + 1, 1, + + 1, 1, + 1, 0, + 1, 1, + + 1, 1, + 1, 0, + 1, 1, + + 1, 1, + 1, 0, + 1, 1, + }; + } + + std::vector actnum_2x3x4_exclude_diagonals() + { + return { + 0, 1, + 1, 1, + 0, 1, + + 1, 1, + 1, 0, + 1, 1, + + 0, 1, + 1, 1, + 0, 1, + + 1, 1, + 1, 0, + 1, 1, + }; + } +} // Namespace Anonymous + +// K = 0 +// +------+------+ +// | 16 | 20 | +// +------+------+ +// | 8 | 12 | +// +------+------+ +// | 0 | 4 | +// +------+------+ +// +// K = 1 +// +------+------+ +// | 17 | 21 | +// +------+------+ +// | 9 | 13 | +// +------+------+ +// | 1 | 5 | +// +------+------+ +// +// K = 2 +// +------+------+ +// | 18 | 22 | +// +------+------+ +// | 10 | 14 | +// +------+------+ +// | 2 | 6 | +// +------+------+ +// +// K = 3 +// +------+------+ +// | 19 | 23 | +// +------+------+ +// | 11 | 15 | +// +------+------+ +// | 3 | 7 | +// +------+------+ +BOOST_AUTO_TEST_CASE(Cube_2x3x4_Full) +{ + const auto grid = Opm::EclipseGrid {{2, 3, 4}, coord_2x3x4(), zcorn_2x3x4() }; + + const auto map = Opm::buildColumnarActiveIndexMappingTables(grid); + + // K = 0 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 0), 0); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 1), 4); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 2), 8); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 3), 12); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 4), 16); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 5), 20); + + // K = 1 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 6), 1); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 7), 5); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 8), 9); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 9), 13); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(10), 17); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(11), 21); + + // K = 2 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(12), 2); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(13), 6); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(14), 10); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(15), 14); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(16), 18); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(17), 22); + + // K = 3 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(18), 3); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(19), 7); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(20), 11); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(21), 15); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(22), 19); + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(23), 23); +} + +// K = 0 +// +------+------+ +// | 12 | 16 | +// +------+------+ +// | :::: | 10 | +// +------+------+ +// | 0 | 4 | +// +------+------+ +// +// K = 1 +// +------+------+ +// | 13 | 17 | +// +------+------+ +// | 8 | :::: | +// +------+------+ +// | 1 | 5 | +// +------+------+ +// +// K = 2 +// +------+------+ +// | 14 | 18 | +// +------+------+ +// | :::: | 11 | +// +------+------+ +// | 2 | 6 | +// +------+------+ +// +// K = 3 +// +------+------+ +// | 15 | 19 | +// +------+------+ +// | 9 | :::: | +// +------+------+ +// | 3 | 7 | +// +------+------+ +BOOST_AUTO_TEST_CASE(Cube_2x3x4_Exclude_Alternate_Centre) +{ + const auto actnum = actnum_2x3x4_exclude_alternate_centre(); + const auto grid = Opm::EclipseGrid {{2, 3, 4}, coord_2x3x4(), zcorn_2x3x4(), actnum.data() }; + + const auto map = Opm::buildColumnarActiveIndexMappingTables(grid); + + // K = 0 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 0), 0); // (0,0,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 1), 4); // (1,0,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 2), 10); // (1,1,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 3), 12); // (0,2,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 4), 16); // (1,2,0) + + // K = 1 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 5), 1); // (0,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 6), 5); // (1,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 7), 8); // (0,1,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 8), 13); // (0,2,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 9), 17); // (1,2,1) + + // K = 2 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(10), 2); // (0,0,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(11), 6); // (1,0,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(12), 11); // (1,1,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(13), 14); // (0,2,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(14), 18); // (1,2,2) + + // K = 3 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(15), 3); // (0,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(16), 7); // (1,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(17), 9); // (0,1,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(18), 15); // (0,2,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(19), 19); // (1,2,3) +} + +// K = 0 +// +------+------+ +// | 12 | 16 | +// +------+------+ +// | 8 | :::: | +// +------+------+ +// | 0 | 4 | +// +------+------+ +// +// K = 1 +// +------+------+ +// | 13 | 17 | +// +------+------+ +// | 9 | :::: | +// +------+------+ +// | 1 | 5 | +// +------+------+ +// +// K = 2 +// +------+------+ +// | 14 | 18 | +// +------+------+ +// | 10 | :::: | +// +------+------+ +// | 2 | 6 | +// +------+------+ +// +// K = 3 +// +------+------+ +// | 15 | 19 | +// +------+------+ +// | 11 | :::: | +// +------+------+ +// | 3 | 7 | +// +------+------+ +BOOST_AUTO_TEST_CASE(Cube_2x3x4_Exclude_Centre_Column) +{ + const auto actnum = actnum_2x3x4_exclude_centre_column(); + const auto grid = Opm::EclipseGrid {{2, 3, 4}, coord_2x3x4(), zcorn_2x3x4(), actnum.data() }; + + const auto map = Opm::buildColumnarActiveIndexMappingTables(grid); + + // K = 0 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 0), 0); // (0,0,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 1), 4); // (1,0,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 2), 8); // (0,1,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 3), 12); // (0,2,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 4), 16); // (1,2,0) + + // K = 1 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 5), 1); // (0,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 6), 5); // (1,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 7), 9); // (0,1,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 8), 13); // (0,2,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 9), 17); // (1,2,1) + + // K = 2 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(10), 2); // (0,0,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(11), 6); // (1,0,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(12), 10); // (0,1,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(13), 14); // (0,2,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(14), 18); // (1,2,2) + + // K = 3 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(15), 3); // (0,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(16), 7); // (1,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(17), 11); // (0,1,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(18), 15); // (0,2,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(19), 19); // (1,2,3) +} + +// K = 0 +// +------+------+ +// | :::: | 14 | +// +------+------+ +// | 6 | 10 | +// +------+------+ +// | :::: | 2 | +// +------+------+ +// +// K = 1 +// +------+------+ +// | 12 | 15 | +// +------+------+ +// | 7 | :::: | +// +------+------+ +// | 0 | 3 | +// +------+------+ +// +// K = 2 +// +------+------+ +// | :::: | 16 | +// +------+------+ +// | 8 | 11 | +// +------+------+ +// | :::: | 4 | +// +------+------+ +// +// K = 3 +// +------+------+ +// | 13 | 17 | +// +------+------+ +// | 9 | :::: | +// +------+------+ +// | 1 | 5 | +// +------+------+ +BOOST_AUTO_TEST_CASE(Cube_2x3x4_Exclude_Diagonals) +{ + const auto actnum = actnum_2x3x4_exclude_diagonals(); + const auto grid = Opm::EclipseGrid {{2, 3, 4}, coord_2x3x4(), zcorn_2x3x4(), actnum.data() }; + + const auto map = Opm::buildColumnarActiveIndexMappingTables(grid); + + // K = 0 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 0), 2); // (0,1,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 1), 6); // (1,0,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 2), 10); // (1,1,0) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 3), 14); // (1,2,0) + + // K = 1 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 4), 0); // (0,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 5), 3); // (1,0,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 6), 7); // (0,1,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 7), 12); // (0,2,1) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 8), 15); // (1,2,1) + + // K = 2 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex( 9), 4); // (1,0,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(10), 8); // (0,1,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(11), 11); // (1,1,2) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(12), 16); // (1,1,2) + + // K = 3 + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(13), 1); // (0,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(14), 5); // (1,0,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(15), 9); // (0,1,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(16), 13); // (0,2,3) + BOOST_CHECK_EQUAL(map.getColumnarActiveIndex(17), 17); // (1,2,3) +} + +BOOST_AUTO_TEST_SUITE_END() // Grid_Based_NY_Larger_Than_NX