diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index f5d7508c..26930e03 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -631,7 +631,7 @@ main(int argc, char** argv) } TransportSolver tsolver(model); // Column-based gravity segregation solver. - typedef std::map > ColMap; + typedef std::pair, std::vector > > ColMap; ColMap columns; if (use_column_solver) { Opm::extractColumn(*grid->c_grid(), columns); diff --git a/opm/core/ColumnExtract.hpp b/opm/core/ColumnExtract.hpp index 1c80a66d..5a4cdc76 100644 --- a/opm/core/ColumnExtract.hpp +++ b/opm/core/ColumnExtract.hpp @@ -32,12 +32,15 @@ struct ExtractColumnCompare /// Extract each column of the grid. /// \note Assumes the pillars of the grid are all vertically aligned. /// \param grid The grid from which to extract the columns. -/// \param columns will for each i + cartgrim[0]*j where (i, j) represents a non-empty column, +/// \param columns will for each (i, j) where (i, j) represents a non-empty column, //// contain the cell indices contained in the column -/// centered at (i, j). -void extractColumn( const UnstructuredGrid& grid, std::map >& columns ) +/// centered at (i, j) in the second variable, and i+jN in the first variable. +void extractColumn( const UnstructuredGrid& grid, std::pair, std::vector > >& columns ) { const int* dims = grid.cartdims; + + // Keeps track of column_index ---> index of vector + std::map global_to_local; for (int i = 0; i < grid.number_of_cells; ++i) { // Extract Cartesian coordinates int index = grid.global_cell ? grid.global_cell[i] : i; // If null, assume mapping is identity. @@ -45,12 +48,22 @@ void extractColumn( const UnstructuredGrid& grid, std::map int k_cart = index / dims[0] / dims[1]; int j_cart = (index - k_cart*dims[0]*dims[1])/ dims[0]; - // Finally add it to the map - columns[i_cart + j_cart*dims[0]].push_back(i); + int local_index; + std::map::iterator local_index_iterator = global_to_local.find(i_cart+j_cart*dims[0]); + if( local_index_iterator != global_to_local.end() ) { + local_index = local_index_iterator->second; + } else { + local_index = columns.second.size(); + global_to_local[i_cart+j_cart*dims[0]] = local_index; + columns.second.push_back(std::vector()); + } + + columns.first.push_back(i_cart + j_cart*dims[0]); + columns.second[local_index].push_back(i); } - for (std::map >::iterator it = columns.begin(); it != columns.end(); ++it) { - std::sort(it->second.begin(), it->second.end(), ExtractColumnCompare(grid)); + for(int i = 0; i < columns.second.size(); ++i) { + std::sort(columns.second[i].begin(), columns.second[i].end(), ExtractColumnCompare(grid)); } } diff --git a/opm/core/transport/GravityColumnSolver.hpp b/opm/core/transport/GravityColumnSolver.hpp index be727ae4..fd1fe0b9 100644 --- a/opm/core/transport/GravityColumnSolver.hpp +++ b/opm/core/transport/GravityColumnSolver.hpp @@ -40,12 +40,12 @@ namespace Opm const double tol, const int maxit); - /// \param[in] columns for each column (with logical cartesian indices as key), + /// \param[in] columns for each column columns.second /// contains the cells on which to solve the segregation /// problem. For each column, its cells must be in a single /// vertical column, and ordered /// (direction doesn't matter). - void solve(const std::map >& columns, + void solve(const std::pair, std::vector > >& columns, const double dt, std::vector& s); private: diff --git a/opm/core/transport/GravityColumnSolver_impl.hpp b/opm/core/transport/GravityColumnSolver_impl.hpp index 0bb32863..a171fd22 100644 --- a/opm/core/transport/GravityColumnSolver_impl.hpp +++ b/opm/core/transport/GravityColumnSolver_impl.hpp @@ -20,7 +20,7 @@ #include #include #include - +#include namespace Opm { @@ -65,14 +65,14 @@ namespace Opm } // anon namespace - - /// \param[in] columns for each column (with logical cartesian indices as key), + + /// \param[in] columns for each column columns.second, /// contains the cells on which to solve the segregation /// problem. For each column, its cells must be in a single /// vertical column, and ordered /// (direction doesn't matter). template - void GravityColumnSolver::solve(const std::map >& columns, + void GravityColumnSolver::solve(const std::pair, std::vector > >& columns, const double dt, std::vector& s) { @@ -84,15 +84,22 @@ namespace Opm int iter = 0; double max_delta = 1e100; + while (iter < maxit_) { model_.initIteration(state, grid_, sys); - std::map >::const_iterator it; - for (it = columns.begin(); it != columns.end(); ++it) { - solveSingleColumn(it->second, dt, s, increment); - } + // std::map >::const_iterator it; + //for (it = columns.begin(); it != columns.end(); ++it) { + int size = columns.second.size(); + + #pragma omp parallel for schedule(dynamic) + for(int i = 0; i < size; ++i) { + solveSingleColumn(columns.second[i], dt, s, increment); + } + for (int cell = 0; cell < grid_.number_of_cells; ++cell) { sys.vector().writableSolution()[cell] += increment[cell]; } + const double maxelem = *std::max_element(increment.begin(), increment.end()); const double minelem = *std::min_element(increment.begin(), increment.end()); max_delta = std::max(maxelem, -minelem); diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 7d9aadac..6c9a5ffd 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -599,7 +599,7 @@ namespace Opm - void TransportModelTwophase::solveGravity(const std::map >& columns, + void TransportModelTwophase::solveGravity(const std::pair, std::vector > >& columns, const double* porevolume, const double dt, std::vector& saturation) @@ -626,10 +626,10 @@ namespace Opm saturation_ = &saturation[0]; // Solve on all columns. - std::map >::const_iterator it; - for (it = columns.begin(); it != columns.end(); ++it) { + + for (int i = 0; i < columns.second.size(); i++) { // std::cout << "==== new column" << std::endl; - solveGravityColumn(it->second); + solveGravityColumn(columns.second[i]); } } diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 8c516898..957a58f2 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -53,7 +53,7 @@ namespace Opm const int pos, const double* gravflux); void solveGravityColumn(const std::vector& cells); - void solveGravity(const std::map >& columns, + void solveGravity(const std::pair, std::vector > >& columns, const double* porevolume, const double dt, std::vector& saturation); diff --git a/tests/test_column_extract.cpp b/tests/test_column_extract.cpp index 2aeb241d..db007c25 100644 --- a/tests/test_column_extract.cpp +++ b/tests/test_column_extract.cpp @@ -23,15 +23,16 @@ BOOST_AUTO_TEST_CASE(SingleColumnTest) grid->global_cell[i] = i; } - std::map > columns; + std::pair, std::vector > > columns; extractColumn(*manager.c_grid(), columns); std::vector correct_answer; + for(int i = 0; i < 10; i++) { correct_answer.push_back(i); } BOOST_CHECK_EQUAL_COLLECTIONS(correct_answer.begin(), correct_answer.end(), - columns[0].begin(), columns[0].end()); + columns.second[0].begin(), columns.second[0].end()); } @@ -49,7 +50,7 @@ BOOST_AUTO_TEST_CASE(FourByFourColumnTest) grid->global_cell[i] = i; } - std::map > columns; + std::pair, std::vector > > columns; extractColumn(*manager.c_grid(), columns); std::vector > correct_answer; @@ -62,7 +63,7 @@ BOOST_AUTO_TEST_CASE(FourByFourColumnTest) for(int i = 0; i < size_x * size_y; i++) { BOOST_CHECK_EQUAL_COLLECTIONS(correct_answer[i].begin(), correct_answer[i].end(), - columns[i].begin(), columns[i].end()); + columns.second[i].begin(), columns.second[i].end()); } }