Make initGravity() to private, call from constructor.

Also modify interface of solveGravity() to be minimal,
construcing columns for segregations solver at construction time.
This commit is contained in:
Atgeirr Flø Rasmussen
2013-03-15 13:53:37 +01:00
parent da985748d0
commit 17ac731a0e
2 changed files with 34 additions and 24 deletions

View File

@@ -21,6 +21,7 @@
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/grid.h>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/utility/RootFinders.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
@@ -42,6 +43,7 @@ namespace Opm
TransportSolverTwophaseReorder::TransportSolverTwophaseReorder(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double* gravity,
const double tol,
const int maxit)
: grid_(grid),
@@ -74,6 +76,10 @@ namespace Opm
cells[i] = i;
}
props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]);
if (gravity) {
initGravity(gravity);
initColumns();
}
}
@@ -552,9 +558,16 @@ namespace Opm
void TransportSolverTwophaseReorder::initColumns()
{
extractColumn(grid_, columns_);
}
void TransportSolverTwophaseReorder::solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux)
const int pos,
const double* gravflux)
{
const int cell = cells[pos];
GravityResidual res(*this, cells, pos, gravflux);
@@ -622,10 +635,9 @@ namespace Opm
void TransportSolverTwophaseReorder::solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation)
void TransportSolverTwophaseReorder::solveGravity(const double* porevolume,
const double dt,
TwophaseState& state)
{
// Initialize mobilities.
const int nc = grid_.number_of_cells;
@@ -634,7 +646,7 @@ namespace Opm
cells[c] = c;
}
mob_.resize(2*nc);
props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0);
props_.relperm(cells.size(), &state.saturation()[0], &cells[0], &mob_[0], 0);
const double* mu = props_.viscosity();
for (int c = 0; c < nc; ++c) {
mob_[2*c] /= mu[0];
@@ -644,18 +656,18 @@ namespace Opm
// Set up other variables.
porevolume_ = porevolume;
dt_ = dt;
toWaterSat(saturation, saturation_);
toWaterSat(state.saturation(), saturation_);
// Solve on all columns.
int num_iters = 0;
for (std::vector<std::vector<int> >::size_type i = 0; i < columns.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[i]);
num_iters += solveGravityColumn(columns_[i]);
}
std::cout << "Gauss-Seidel column solver average iterations: "
<< double(num_iters)/double(columns.size()) << std::endl;
<< double(num_iters)/double(columns_.size()) << std::endl;
toBothSat(saturation_, saturation);
toBothSat(saturation_, state.saturation());
}
} // namespace Opm

View File

@@ -43,6 +43,7 @@ namespace Opm
/// \param[in] maxit Maximum number of non-linear iterations used.
TransportSolverTwophaseReorder(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double* gravity,
const double tol,
const int maxit);
@@ -60,29 +61,25 @@ namespace Opm
const double dt,
TwophaseState& state);
/// Initialise quantities needed by gravity solver.
/// \param[in] grav gravity vector
void initGravity(const double* grav);
/// Solve for gravity segregation.
/// This uses a column-wise nonlinear Gauss-Seidel approach.
/// It assumes that the input columns contain cells in a single
/// vertical stack, that do not interact with other columns (for
/// gravity segregation.
/// \param[in] columns Vector of cell-columns.
/// It assumes that the grid can be divided into vertical columns
/// that do not interact with each other (for gravity segregation).
/// \param[in] porevolume Array of pore volumes.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
/// \param[in, out] state Reservoir state. Calling solveGravity() will read state.faceflux() and
/// read and write state.saturation().
void solveGravity(const double* porevolume,
const double dt,
std::vector<double>& saturation);
TwophaseState& state);
//// Return the number of iterations used by the reordering solver.
//// \return vector of iteration per cell
const std::vector<int>& getReorderIterations() const;
private:
void initGravity(const double* grav);
void initColumns();
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
@@ -111,6 +108,7 @@ namespace Opm
std::vector<double> gravflux_;
std::vector<double> mob_;
std::vector<double> s0_;
std::vector<std::vector<int> > columns_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;