OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make SimulatorIncompTwophase flexible w.r.t. transport solver.

Browse Source 
It can now use the reordering or the full Newton-Raphson solvers.
This commit is contained in:
Atgeirr Flø Rasmussen
2013-03-15 14:16:59 +01:00
parent b8ed52d89e
commit 25bacfd231
1 changed files with 42 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
opm/simulators/SimulatorIncompTwophase.cpp
View File

@@ -42,10 +42,9 @@
#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>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/transport/TransportSolverTwophaseImplicit.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
@@ -87,6 +86,7 @@ namespace Opm
int max_well_control_iterations_;
// Parameters for transport solver.
int num_transport_substeps_;
bool use_reorder_;
bool use_segregation_split_;
// Observed objects.
const UnstructuredGrid& grid_;
@@ -98,11 +98,7 @@ namespace Opm
const FlowBoundaryConditions* bcs_;
// Solvers
IncompTpfa psolver_;
// this should maybe be packed in a separate file
boost::scoped_ptr<TransportSolverTwophaseInterface> tsolver_;
//ImpliciteTwoPhaseTransportSolver tsolver_;
// Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_;
// Misc. data
std::vector<int> allcells_;
};
@@ -320,7 +316,9 @@ namespace Opm
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
: grid_(grid),
: use_reorder_(param.getDefault("use_reorder", true)),
use_segregation_split_(param.getDefault("use_segregation_split", false)),
grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
@@ -333,23 +331,24 @@ namespace Opm
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs)
{
const bool use_reorder = param.getDefault("use_reorder", true);
if (use_reorder) {
/*
tsolver_.reset(new Opm::TransportModelTwoPhase(grid, props, rock_comp_props, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs));
*/
// Initialize transport solver.
if (use_reorder_) {
tsolver_.reset(new Opm::TransportSolverTwophaseReorder(grid,
props,
use_segregation_split_ ? gravity : NULL,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)));
} else {
if (rock_comp_props->isActive()) {
THROW("The implicit pressure solver cannot handle rock compressibility.");
}
if (use_segregation_split_) {
THROW("The implicit pressure solver is not set up to use segregation splitting.");
}
std::vector<double> porevol;
computePorevolume(grid, props.porosity(), porevol);
tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(
grid,
tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(grid,
props,
porevol,
gravity,
@@ -379,7 +378,7 @@ namespace Opm
// Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
// Misc init.
const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells);
@@ -445,6 +444,12 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
if (use_reorder_) {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
}
SimulatorReport sreport;
@@ -538,10 +543,7 @@ namespace Opm
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_->solve(&transport_src[0],
&porevol[0],
stepsize,
state);
tsolver_->solve(&transport_src[0], &porevol[0], stepsize, state);
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
@@ -551,18 +553,20 @@ namespace Opm
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (use_reorder_ && use_segregation_split_) {
// Again, unfortunate but safe use of dynamic_cast.
// Possible solution: refactor gravity solver to its own class.
dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
.solveGravity(&initial_porevol[0], stepsize, state);
}
watercut.push(timer.currentTime() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells_) {
wellreport.push(props_, *wells_, state.saturation(),
timer.currentTime() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
@@ -590,6 +594,12 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
if (use_reorder_) {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
outputWaterCut(watercut, output_dir_);
if (wells_) {
outputWellReport(wellreport, output_dir_);