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

Make initGravity() to private, call from constructor.

Browse Source 
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 aadae49b41
commit b304105b4f
6 changed files with 62 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
examples/spu_2p.cpp
View File

@ -455,10 +455,9 @@ main(int argc, char** argv)
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportSolverTwophaseReorder reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_gauss_seidel_gravity) {
reorder_model.initGravity(grav);
}
const double* transport_grav = use_gauss_seidel_gravity ? grav : NULL;
Opm::TransportSolverTwophaseReorder reorder_model(*grid->c_grid(), *props, transport_grav, nl_tolerance, nl_maxiter);
// Non-reordering solver.
TransportModel model (fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
if (use_gravity) {
@ -633,7 +632,7 @@ main(int argc, char** argv)
if (use_segregation_split) {
if (use_column_solver) {
if (use_gauss_seidel_gravity) {
reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
reorder_model.solveGravity(&porevol[0], stepsize, state);
} else {
colsolver.solve(columns, stepsize, state.saturation());
}

49
opm/core/simulator/SimulatorIncompTwophaseReorder.cpp
View File

@ -42,7 +42,6 @@
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
@ -109,14 +108,14 @@ namespace Opm
SimulatorIncompTwophaseReorder::SimulatorIncompTwophaseReorder(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity));
}
@ -126,8 +125,8 @@ namespace Opm
SimulatorReport SimulatorIncompTwophaseReorder::run(SimulatorTimer& timer,
TwophaseState& state,
WellState& well_state)
TwophaseState& state,
WellState& well_state)
{
return pimpl_->run(timer, state, well_state);
}
@ -310,15 +309,16 @@ namespace Opm
SimulatorIncompTwophaseReorder::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
: grid_(grid),
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
: use_segregation_split_(param.getDefault("use_segregation_split", false)),
grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
@ -330,7 +330,7 @@ namespace Opm
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs),
tsolver_(grid, props,
tsolver_(grid, props, use_segregation_split_ ? gravity : NULL,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30))
{
@ -356,11 +356,6 @@ namespace Opm
// Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
// Misc init.
const int num_cells = grid.number_of_cells;
@ -533,7 +528,7 @@ namespace Opm
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (use_segregation_split_) {
tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
tsolver_.solveGravity(&initial_porevol[0], stepsize, state);
}
watercut.push(timer.currentTime() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),

36
opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp
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

22
opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp
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_;

2
tutorials/tutorial3.cpp
View File

@ -232,7 +232,7 @@ int main ()
/// \internal [transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, tolerance, max_iterations);
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, NULL, tolerance, max_iterations);
/// \internal [transport solver]
/// \endinternal

2
tutorials/tutorial4.cpp
View File

@ -169,7 +169,7 @@ int main ()
/// \internal[transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, tolerance, max_iterations);
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, NULL, tolerance, max_iterations);
/// \internal[transport solver]
/// \endinternal