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

Simplified data structure for extracted columns.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-05-11 09:37:12 +02:00
parent c521c0a0f6
commit 891696e333
7 changed files with 27 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
examples/spu_2p.cpp
View File

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

22
opm/core/ColumnExtract.hpp
View File

@ -35,36 +35,34 @@ struct ExtractColumnCompare
/// \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) in the second variable, and i+jN in the first variable.
inline void extractColumn( const UnstructuredGrid& grid, std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns )
inline void extractColumn( const UnstructuredGrid& grid, 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) {
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
// Extract Cartesian coordinates
int index = grid.global_cell ? grid.global_cell[i] : i; // If null, assume mapping is identity.
int index = grid.global_cell ? grid.global_cell[cell] : cell; // If null, assume mapping is identity.
int i_cart = index % dims[0];
int k_cart = index / dims[0] / dims[1];
int j_cart = (index - k_cart*dims[0]*dims[1])/ dims[0];
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() ) {
if (local_index_iterator != global_to_local.end()) {
local_index = local_index_iterator->second;
} else {
local_index = columns.second.size();
local_index = columns.size();
global_to_local[i_cart+j_cart*dims[0]] = local_index;
columns.second.push_back(std::vector<int>());
columns.push_back(std::vector<int>());
}
columns.first.push_back(i_cart + j_cart*dims[0]);
columns.second[local_index].push_back(i);
columns[local_index].push_back(cell);
}
int num_cols = columns.second.size();
for(int i = 0; i < num_cols; ++i) {
std::sort(columns.second[i].begin(), columns.second[i].end(), ExtractColumnCompare(grid));
int num_cols = columns.size();
for (int col = 0; col < num_cols; ++col) {
std::sort(columns[col].begin(), columns[col].end(), ExtractColumnCompare(grid));
}
}

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

@ -40,16 +40,16 @@ namespace Opm
const double tol,
const int maxit);
/// \param[in] columns for each column columns.second
/// \param[in] columns for each column col, columns[col]
/// contains the cells on which to solve the segregation
/// problem. For each column, its cells must be in a single
/// vertical column, and ordered
/// vertical column, connected and ordered
/// (direction doesn't matter).
void solve(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
void solve(const std::vector<std::vector<int> >& columns,
const double dt,
std::vector<double>& s);
private:
void solveSingleColumn(const std::vector<int>& column_cells,
void solveSingleColumn(const std::vector<int>& columns,
const double dt,
std::vector<double>& s,
std::vector<double>& sol_vec);

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

@ -72,7 +72,7 @@ namespace Opm
/// vertical column, and ordered
/// (direction doesn't matter).
template <class Model>
void GravityColumnSolver<Model>::solve(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
void GravityColumnSolver<Model>::solve(const std::vector<std::vector<int> >& columns,
const double dt,
std::vector<double>& s)
{
@ -89,11 +89,11 @@ namespace Opm
model_.initIteration(state, grid_, sys);
// std::map<int, std::vector<int> >::const_iterator it;
//for (it = columns.begin(); it != columns.end(); ++it) {
int size = columns.second.size();
int size = columns.size();
#pragma omp parallel for schedule(dynamic)
for(int i = 0; i < size; ++i) {
solveSingleColumn(columns.second[i], dt, s, increment);
solveSingleColumn(columns[i], dt, s, increment);
}
for (int cell = 0; cell < grid_.number_of_cells; ++cell) {

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

@ -602,7 +602,7 @@ namespace Opm
void TransportModelTwophase::solveGravity(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
void TransportModelTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation)
@ -630,12 +630,12 @@ namespace Opm
// Solve on all columns.
int num_iters = 0;
for (std::vector<std::vector<int> >::size_type i = 0; i < columns.second.size(); i++) {
for (std::vector<std::vector<int> >::size_type i = 0; i < columns.size(); i++) {
// std::cout << "==== new column" << std::endl;
num_iters += solveGravityColumn(columns.second[i]);
num_iters += solveGravityColumn(columns[i]);
}
std::cout << "Gauss-Seidel column solver average iterations: "
<< double(num_iters)/double(columns.second.size()) << std::endl;
<< double(num_iters)/double(columns.size()) << std::endl;
}
} // namespace Opm

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

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

8
tests/test_column_extract.cpp
View File

@ -23,7 +23,7 @@ BOOST_AUTO_TEST_CASE(SingleColumnTest)
grid->global_cell[i] = i;
}
std::pair<std::vector<int>, std::vector<std::vector<int> > > columns;
std::vector<std::vector<int> > columns;
extractColumn(*manager.c_grid(), columns);
std::vector<int> correct_answer;
@ -32,7 +32,7 @@ BOOST_AUTO_TEST_CASE(SingleColumnTest)
correct_answer.push_back(i);
}
BOOST_CHECK_EQUAL_COLLECTIONS(correct_answer.begin(), correct_answer.end(),
columns.second[0].begin(), columns.second[0].end());
columns[0].begin(), columns[0].end());
}
@ -50,7 +50,7 @@ BOOST_AUTO_TEST_CASE(FourByFourColumnTest)
grid->global_cell[i] = i;
}
std::pair<std::vector<int>, std::vector<std::vector<int> > > columns;
std::vector<std::vector<int> > columns;
extractColumn(*manager.c_grid(), columns);
std::vector<std::vector<int> > correct_answer;
@ -63,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.second[i].begin(), columns.second[i].end());
columns[i].begin(), columns[i].end());
}
}