OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added opm-parallelization to column solver for gravity. Did a (minor) interface change on extractColumn

Browse Source
This commit is contained in:
Kjetil Olsen Lye 2012-03-23 12:32:25 +01:00
parent 355eb053e8
commit dd2d474643
7 changed files with 48 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/spu_2p.cpp
View File

@ -631,7 +631,7 @@ main(int argc, char** argv)
}
TransportSolver tsolver(model);
// Column-based gravity segregation solver.
typedef std::map<int, std::vector<int> > ColMap;
typedef std::pair<std::vector<int>, std::vector<std::vector<int> > > ColMap;
ColMap columns;
if (use_column_solver) {
Opm::extractColumn(*grid->c_grid(), columns);

27
opm/core/ColumnExtract.hpp
View File

@ -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<int, std::vector<int> >& 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<int>, std::vector<std::vector<int> > >& columns )
{
const int* dims = grid.cartdims;
// Keeps track of column_index ---> index of vector
std::map<int, int> 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, std::vector<int>
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<int, int>::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<int>());
}
columns.first.push_back(i_cart + j_cart*dims[0]);
columns.second[local_index].push_back(i);
}
for (std::map<int, std::vector<int> >::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));
}
}

4
opm/core/transport/GravityColumnSolver.hpp
View File

@ -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<int, std::vector<int> >& columns,
void solve(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
const double dt,
std::vector<double>& s);
private:

23
opm/core/transport/GravityColumnSolver_impl.hpp
View File

@ -20,7 +20,7 @@
#include <opm/core/transport/GravityColumnSolver.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <sys/time.h>
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 <class Model>
void GravityColumnSolver<Model>::solve(const std::map<int, std::vector<int> >& columns,
void GravityColumnSolver<Model>::solve(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
const double dt,
std::vector<double>& s)
{
@ -84,15 +84,22 @@ namespace Opm
int iter = 0;
double max_delta = 1e100;
while (iter < maxit_) {
model_.initIteration(state, grid_, sys);
std::map<int, std::vector<int> >::const_iterator it;
for (it = columns.begin(); it != columns.end(); ++it) {
solveSingleColumn(it->second, dt, s, increment);
}
// std::map<int, std::vector<int> >::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);

8
opm/core/transport/reorder/TransportModelTwophase.cpp
View File

@ -599,7 +599,7 @@ namespace Opm
void TransportModelTwophase::solveGravity(const std::map<int, std::vector<int> >& columns,
void TransportModelTwophase::solveGravity(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation)
@ -626,10 +626,10 @@ namespace Opm
saturation_ = &saturation[0];
// Solve on all columns.
std::map<int, std::vector<int> >::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]);
}
}

2
opm/core/transport/reorder/TransportModelTwophase.hpp
View File

@ -53,7 +53,7 @@ namespace Opm
const int pos,
const double* gravflux);
void solveGravityColumn(const std::vector<int>& cells);
void solveGravity(const std::map<int, std::vector<int> >& columns,
void solveGravity(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation);

9
tests/test_column_extract.cpp
View File

@ -23,15 +23,16 @@ BOOST_AUTO_TEST_CASE(SingleColumnTest)
grid->global_cell[i] = i;
}
std::map<int, std::vector<int> > columns;
std::pair<std::vector<int>, std::vector<std::vector<int> > > columns;
extractColumn(*manager.c_grid(), columns);
std::vector<int> 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<int, std::vector<int> > columns;
std::pair<std::vector<int>, std::vector<std::vector<int> > > columns;
extractColumn(*manager.c_grid(), columns);
std::vector<std::vector<int> > 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());
}
}