Changes needed to for makeing a simulator using ImplicitTransport. Several changes in names to highlight what is reorder simulator classes

This commit is contained in:
Halvor Møll Nilsen 2012-11-16 13:38:03 +01:00
parent c090eb618b
commit 523ca56ef1

View File

@ -45,8 +45,9 @@
#include <opm/core/utility/ColumnExtract.hpp> #include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/TwophaseState.hpp> #include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
//#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp> #include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/transport/ImpliciteTwoPhaseTransportSolver.hpp>
#include <boost/filesystem.hpp> #include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp> #include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp> #include <boost/lexical_cast.hpp>
@ -99,28 +100,8 @@ namespace Opm
// Solvers // Solvers
IncompTpfa psolver_; IncompTpfa psolver_;
// this should maybe be packed in a separate file // this should maybe be packed in a separate file
typedef Opm::SimpleFluid2pWrappingProps TwophaseFluid; boost::scoped_ptr<TwoPhaseTransportSolver> tsolver_;
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel; //ImpliciteTwoPhaseTransportSolver tsolver_;
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return AccumulationNorm <Vector, MaxAbs>::norm(v);
}
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > ImpliciteTwoPhaseTransportSolver;
ImpliciteTwoPhaseTransportSolver tsolver_;
// Needed by column-based gravity segregation solver. // Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_; std::vector< std::vector<int> > columns_;
// Misc. data // Misc. data
@ -339,7 +320,6 @@ namespace Opm
const std::vector<double>& src, const std::vector<double>& src,
const FlowBoundaryConditions* bcs, const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolverp, LinearSolverInterface& linsolverp,
LinearSolverInterface& linsolvert,
const double* gravity) const double* gravity)
: grid_(grid), : grid_(grid),
props_(props), props_(props),
@ -352,23 +332,41 @@ namespace Opm
param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10), param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs), gravity, wells_manager.c_wells(), src, bcs)
{ {
const bool use_reorder = param.getDefault("use_reorder", true); const bool use_reorder = param.getDefault("use_reorder", true);
if(use_reorder){ if(use_reorder){
tsolver_ = new TransportSolverTwoPhaseReorder(grid, props, rock_comp_props, linsolver, /*
param.getDefault("nl_pressure_residual_tolerance", 0.0), tsolver_.reset(new Opm::TransportModelTwoPhase(grid, props, rock_comp_props, linsolver,
param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_maxiter", 10), param.getDefault("nl_pressure_change_tolerance", 1.0),
gravity, wells_manager.c_wells(), src, bcs); param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs));
*/
}else{ }else{
//Opm::ImplicitTransportDetails::NRReport rpt; //Opm::ImplicitTransportDetails::NRReport rpt;
Opm::ImplicitTransportDetails::NRControl ctrl; Opm::ImplicitTransportDetails::NRControl ctrl;
ctrl.max_it = param.getDefault("max_it", 20); ctrl.max_it = param.getDefault("max_it", 20);
ctrl.verbosity = param.getDefault("verbosity", 0); ctrl.verbosity = param.getDefault("verbosity", 0);
ctrl.max_it_ls = param.getDefault("max_it_ls", 5); ctrl.max_it_ls = param.getDefault("max_it_ls", 5);
tsolver_ = new ImpliciteTransportSolverTwoPhase(grid, props, const bool guess_old_solution = param.getDefault("guess_old_solution", false);
ctrl); Opm::SimpleFluid2pWrappingProps fluid(props);
std::vector<double> porevol;
//if (rock_comp->isActive()) {
// computePorevolume(grid, props->porosity(), rock_comp, state.pressure(), porevol);
//} else {
computePorevolume(grid, props.porosity(), porevol);
//}
SinglePointUpwindTwoPhase<Opm::SimpleFluid2pWrappingProps>
model(fluid, grid, porevol, gravity, guess_old_solution);
tsolver_.reset(new Opm::ImpliciteTwoPhaseTransportSolver(
wells_manager,
*rock_comp_props,
ctrl,
model,
grid,
props,
param));
} }
@ -394,12 +392,7 @@ namespace Opm
// Transport related init. // Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false); use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
// Misc init. // Misc init.
const int num_cells = grid.number_of_cells; const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells); allcells_.resize(num_cells);
@ -465,9 +458,6 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
} }
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputVectorMatlab(std::string("reorder_it"),
tsolver_.getReorderIterations(),
timer.currentStepNum(), output_dir_);
} }
SimulatorReport sreport; SimulatorReport sreport;
@ -563,7 +553,12 @@ namespace Opm
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
//tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], //tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0],
// stepsize, state.saturation()); // stepsize, state.saturation());
tsolver_.solve(*grid_->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt); tsolver_->solve(&transport_src[0],
&porevol[0],
stepsize,
state,
well_state);
double substep_injected[2] = { 0.0 }; double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 }; double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize,
@ -611,9 +606,6 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
} }
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputVectorMatlab(std::string("reorder_it"),
tsolver_.getReorderIterations(),
timer.currentStepNum(), output_dir_);
outputWaterCut(watercut, output_dir_); outputWaterCut(watercut, output_dir_);
if (wells_) { if (wells_) {
outputWellReport(wellreport, output_dir_); outputWellReport(wellreport, output_dir_);