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 aaa10cf17b
commit 2003ea5887
7 changed files with 201 additions and 158 deletions

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_);

View File

@ -27,50 +27,59 @@
*/ */
#include "ImpliciteTwoPhaseTransportSolver.hpp" #include <opm/core/transport/ImpliciteTwoPhaseTransportSolver.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
namespace Opm{
ImpliciteTwoPhaseTransportSolver::ImpliciteTwoPhaseTransportSolver() ImpliciteTwoPhaseTransportSolver::ImpliciteTwoPhaseTransportSolver(
{ const Opm::WellsManager& wells,
ImpliciteTwoPhaseTransportSolver::ImpliciteTwoPhaseTransportSolver(TwoPhaseTransprotModel& model, const Opm::RockCompressibility& rock_comp,
const Opm::IncompPropertiesInterface& param) const ImplicitTransportDetails::NRControl& ctrl,
SinglePointUpwindTwoPhase<Opm::SimpleFluid2pWrappingProps>& model,
const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const parameter::ParameterGroup& param)
: tsolver_(model), : tsolver_(model),
grid_(grid), grid_(grid),
ctrl_(ctrl), ctrl_(ctrl),
props_(props), props_(props),
rock_comp_(rock_comp), rock_comp_(rock_comp),
wells_(wells), wells_(wells)
{ {
tsrc_ = create_transport_source(2, 2); tsrc_ = create_transport_source(2, 2);
//linsolver_(param),
//src_(num_cells, 0.0); //src_(num_cells, 0.0);
} }
void ImpliciteTwoPhaseTransportSolver::solve(const double* darcyflux, void ImpliciteTwoPhaseTransportSolver::solve(const double* porevolume,
const double* porevolume, const double* source,
const double* source, const double dt,
const double dt, TwophaseState& state,
std::vector<double>& saturation) WellState& well_state)
{ {
if (rock_comp->isActive()) { std::vector<double> porevol;
computePorevolume(grid_->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); if (rock_comp_.isActive()) {
computePorevolume(grid_, props_.porosity(), rock_comp_, state.pressure(), porevol);
} }
std::vector<double> src(num_cells, 0.0); std::vector<double> src(grid_.number_of_cells, 0.0);
//Opm::wellsToSrc(*wells->c_wells(), num_cells, src); //Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, Opm::computeTransportSource(grid_, src, state.faceflux(), 1.0,
wells->c_wells(), well_state.perfRates(), src); wells_.c_wells(), well_state.perfRates(), src);
double ssrc[] = { 1.0, 0.0 }; double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 }; double ssink[] = { 0.0, 1.0 };
double zdummy[] = { 0.0, 0.0 }; double zdummy[] = { 0.0, 0.0 };
for (int cell = 0; cell < num_cells; ++cell) { for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
clear_transport_source(tsrc); clear_transport_source(tsrc_);
if (src[cell] > 0.0) { if (src[cell] > 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc); append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc_);
} else if (src[cell] < 0.0) { } else if (src[cell] < 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc_);
} }
} }
// Boundary conditions. // Boundary conditions.
Opm::ImplicitTransportDetails::NRReport rpt; Opm::ImplicitTransportDetails::NRReport rpt;
tsolver.solve(grid_->c_grid(), tsrc, dt, ctrl_, state, linsolvet_, rpt); tsolver_.solve(grid_, tsrc_, dt, ctrl_, state, linsolver_, rpt);
std::cout << rpt; std::cout << rpt;
} }

View File

@ -25,69 +25,103 @@
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <vector>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#ifndef IMPLICITETWOPHASETRANSPORTSOLVER_HPP #ifndef IMPLICITETWOPHASETRANSPORTSOLVER_HPP
#define IMPLICITETWOPHASETRANSPORTSOLVER_HPP #define IMPLICITETWOPHASETRANSPORTSOLVER_HPP
#include <vector>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/transport/TwoPhaseTransportSolver.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/transport/SimpleFluid2pWrappingProps.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/transport/ImplicitTransport.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverFactory.hpp>
// implicite transprot solver #include <opm/core/transport/transport_source.h>
class ImpliciteTwoPhaseTransportSolver : public TwoPhaseTransportSolver #include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
{ #include <opm/core/transport/NormSupport.hpp>
public: #include <opm/core/transport/ImplicitAssembly.hpp>
/// Construct solver. #include <opm/core/transport/ImplicitTransport.hpp>
/// \param[in] grid A 2d or 3d grid. #include <opm/core/transport/JacobianSystem.hpp>
/// \param[in] props Rock and fluid properties. #include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
/// \param[in] tol Tolerance used in the solver. #include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
/// \param[in] maxit Maximum number of non-linear iterations used. #include <boost/scoped_ptr.hpp>
TransportModelTwophase(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
Opm::parameter::ParameterGroup& param);
#include <opm/core/fluid/RockCompressibility.hpp>
/// Solve for saturation at next timestep. #include <opm/core/wells/WellsManager.hpp>
/// \param[in] darcyflux Array of signed face fluxes. #include <opm/core/simulator/WellState.hpp>
/// \param[in] porevolume Array of pore volumes. namespace Opm{
/// \param[in] source Transport source term. // implicite transprot solver
/// \param[in] dt Time step. class ImpliciteTwoPhaseTransportSolver : public TwoPhaseTransportSolver
/// \param[in, out] saturation Phase saturations. {
void solve(const double* darcyflux,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation);
private:
typedef SimpleFluid2pWrappingProps TwophaseFluid;
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public: public:
static double /// Construct solver.
norm(const Vector& v) { /// \param[in] grid A 2d or 3d grid.
return AccumulationNorm <Vector, MaxAbs>::norm(v); /// \param[in] props Rock and fluid properties.
} /// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used.
ImpliciteTwoPhaseTransportSolver(
const Opm::WellsManager& wells,
const Opm::RockCompressibility& rock_comp,
const ImplicitTransportDetails::NRControl& ctrl,
SinglePointUpwindTwoPhase<Opm::SimpleFluid2pWrappingProps>& model,
const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const parameter::ParameterGroup& param);
//ImpliciteTwoPhaseTransportSolver(){
// destroy_transport_source(tsrc_);
//}
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
void solve(const double* porevolume,
const double* source,
const double dt,
Opm::TwophaseState& state,
Opm::WellState& well_state);
private:
typedef Opm::SimpleFluid2pWrappingProps TwophaseFluid;
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
//using namespace ImplicitTransportDefault;
typedef ImplicitTransportDefault::NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef ImplicitTransportDefault::JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return ImplicitTransportDefault::AccumulationNorm <Vector, ImplicitTransportDefault::MaxAbs>::norm(v);
}
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
ImplicitTransportDefault::VectorNegater ,
ImplicitTransportDefault::VectorZero ,
ImplicitTransportDefault::MatrixZero ,
ImplicitTransportDefault::VectorAssign > TransportSolver;
//Opm::LinearSolverFactory
Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolver_;
TransportSolver tsolver_;
const UnstructuredGrid& grid_;
const Opm::ImplicitTransportDetails::NRControl& ctrl_;
const Opm::IncompPropertiesInterface& props_;
const Opm::RockCompressibility& rock_comp_;
const Opm::WellsManager& wells_;
TransportSource* tsrc_;
}; };
}
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
TransportSolver tsolver_,
UnstructuredGrid* grid_,
Opm::ImplicitTransportDetails::NRControl ctrl_,
boost::scoped_ptr<Opm::IncompPropertiesInterface> props_,
boost::scoped_ptr<Opm::RockCompressibility> rock_comp_(rock_comp),
boost::scoped_ptr<Opm::WellsManager> wells_(wells),
};
#endif // IMPLICITETWOPHASETRANSPORTSOLVER_HPP #endif // IMPLICITETWOPHASETRANSPORTSOLVER_HPP

View File

@ -27,17 +27,15 @@
*/ */
#include "SimpleFluid2pWrappingProps.hpp" #include <opm/core/transport/SimpleFluid2pWrappingProps.hpp>
#include <cassert>
SimpleFluid2pWrappingProps::SimpleFluid2pWrappingProps() #include <opm/core/utility/ErrorMacros.hpp>
{ namespace Opm{
}
namespace opm{
SimpleFluid2pWrappingProps::SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props) SimpleFluid2pWrappingProps::SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props)
: props_(props), : props_(props),
smin_(props.numCells()*props.numPhases()), smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases()) smax_(props.numCells()*props.numPhases())
{ {
if (props.numPhases() != 2) { if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases."); THROW("SimpleFluid2pWrapper requires 2 phases.");
@ -50,7 +48,7 @@ namespace opm{
props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]);
} }
double density(int phase) const double SimpleFluid2pWrappingProps::density(int phase) const
{ {
return props_.density()[phase]; return props_.density()[phase];
} }
@ -58,7 +56,7 @@ namespace opm{
template <class Sat, template <class Sat,
class Mob, class Mob,
class DMob> class DMob>
void SimpleFluid2pWrappingProps:;mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const void SimpleFluid2pWrappingProps::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{ {
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity(); const double* mu = props_.viscosity();

View File

@ -29,7 +29,11 @@
#ifndef SIMPLEFLUID2PWRAPPINGPROPS_HPP #ifndef SIMPLEFLUID2PWRAPPINGPROPS_HPP
#define SIMPLEFLUID2PWRAPPINGPROPS_HPP #define SIMPLEFLUID2PWRAPPINGPROPS_HPP
namespace opm{
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <vector>
namespace Opm{
class SimpleFluid2pWrappingProps class SimpleFluid2pWrappingProps
{ {
public: public:
@ -59,7 +63,7 @@ namespace opm{
const Opm::IncompPropertiesInterface& props_; const Opm::IncompPropertiesInterface& props_;
std::vector<double> smin_; std::vector<double> smin_;
std::vector<double> smax_; std::vector<double> smax_;
} };
} }

View File

@ -29,6 +29,4 @@
#include "TwoPhaseTransportSolver.hpp" #include "TwoPhaseTransportSolver.hpp"
TwoPhaseTransportSolver::TwoPhaseTransportSolver()
{
}

View File

@ -29,24 +29,32 @@
#ifndef TWOPHASETRANSPORTSOLVER_HPP #ifndef TWOPHASETRANSPORTSOLVER_HPP
#define TWOPHASETRANSPORTSOLVER_HPP #define TWOPHASETRANSPORTSOLVER_HPP
/// Base class for tranport solvers
class TwoPhaseTransportSolver
{
public:
/// Construct solver.
virtual ~TwoPhaseTransportSolver() {};
/// Solve for saturation at next timestep. #include <opm/core/simulator/TwophaseState.hpp>
/// \param[in] darcyflux Array of signed face fluxes. #include <opm/core/simulator/WellState.hpp>
/// \param[in] porevolume Array of pore volumes. namespace Opm
/// \param[in] source Transport source term. {
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations. /// Base class for tranport solvers
virtual void solve(const double* darcyflux, class TwoPhaseTransportSolver
const double* porevolume, {
const double* source, public:
const double dt, /// Virtual destructor to enable inheritance.
std::vector<double>& saturation) = 0; virtual ~TwoPhaseTransportSolver() {}
};
/// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
virtual void solve(const double* porevolume,
const double* source,
const double dt,
TwophaseState& state,
WellState& wstate) = 0;
};
}
#endif // TWOPHASETRANSPORTSOLVER_HPP #endif // TWOPHASETRANSPORTSOLVER_HPP