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

properly use the SimulatorBase class for SimulatorFullyImplicitBlackoilPolymer and SimulatorFullyImplicitCompressiblePolymer

Browse Source
This commit is contained in:
Andreas Lauser 2015-05-27 23:15:30 +02:00
parent 0f1a7a16d7
commit d27fb2bc45
9 changed files with 409 additions and 264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
examples/flow_polymer.cpp
View File

@ -179,10 +179,10 @@ try
grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
Opm::PolymerBlackoilOutputWriter outputWriter(cGrid,
param,
eclipseState,
pu );
Opm::BlackoilOutputWriter outputWriter(cGrid,
param,
eclipseState,
pu );
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
@ -201,8 +201,8 @@ try
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state.blackoilState());
initBlackoilSurfvol(*grid->c_grid(), *props, state.blackoilState());
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int np = props->numPhases();
@ -215,10 +215,10 @@ try
} else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
state.init(*grid->c_grid(), props->numPhases());
const double grav = param.getDefault("gravity", unit::gravity);
initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState());
initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state);
state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
} else {
initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState());
initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
@ -268,21 +268,21 @@ try
std::vector<double> threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid());
SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid> simulator(param,
*grid->c_grid(),
geology,
*new_props,
polymer_props_ad,
rock_comp->isActive() ? rock_comp.get() : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
polymer,
eclipseState,
outputWriter,
deck,
threshold_pressures);
SimulatorFullyImplicitBlackoilPolymer simulator(param,
*grid->c_grid(),
geology,
*new_props,
polymer_props_ad,
rock_comp->isActive() ? rock_comp.get() : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
polymer,
eclipseState,
outputWriter,
deck,
threshold_pressures);
if (!schedule->initOnly()){
std::cout << "\n\n================ Starting main simulation loop ===============\n"

10
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -167,10 +167,10 @@ try
grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
Opm::PolymerBlackoilOutputWriter outputWriter(cGrid,
param,
eclipseState,
pu );
Opm::BlackoilOutputWriter outputWriter(cGrid,
param,
eclipseState,
pu );
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
@ -224,7 +224,7 @@ try
SimulatorReport fullReport;
// Create and run simulator.
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
SimulatorFullyImplicitCompressiblePolymer<UnstructuredGrid>
SimulatorFullyImplicitCompressiblePolymer
simulator(param,
*grid->c_grid(),
geology,

2
opm/polymer/CompressibleTpfaPolymer.cpp
View File

@ -90,7 +90,7 @@ namespace Opm
{
c_ = &state.concentration();
cmax_ = &state.maxconcentration();
CompressibleTpfa::solve(dt, state.blackoilState(), well_state);
CompressibleTpfa::solve(dt, state, well_state);
}
/// Compute per-solve dynamic properties.

39
opm/polymer/PolymerBlackoilState.hpp
View File

@ -30,7 +30,7 @@ namespace Opm
/// Simulator state for a compressible two-phase simulator with polymer.
/// We use the Blackoil state parameters.
class PolymerBlackoilState
class PolymerBlackoilState : public BlackoilState
{
public:
void init(const UnstructuredGrid& g, int num_phases)
@ -40,53 +40,18 @@ namespace Opm
void init(int number_of_cells, int number_of_faces, int num_phases)
{
state_blackoil_.init(number_of_cells, number_of_faces, num_phases);
BlackoilState::init(number_of_cells, number_of_faces, num_phases);
concentration_.resize(number_of_cells, 0.0);
cmax_.resize(number_of_cells, 0.0);
}
int numPhases() const
{
return state_blackoil_.numPhases();
}
enum ExtremalSat { MinSat = BlackoilState::MinSat, MaxSat = BlackoilState::MaxSat };
void setFirstSat(const std::vector<int>& cells,
const Opm::BlackoilPropertiesInterface& props,
ExtremalSat es)
{
// A better solution for embedding BlackoilState::ExtremalSat could perhaps
// be found, to avoid the cast.
state_blackoil_.setFirstSat(cells, props, static_cast<BlackoilState::ExtremalSat>(es));
}
std::vector<double>& pressure () { return state_blackoil_.pressure(); }
std::vector<double>& temperature () { return state_blackoil_.temperature(); }
std::vector<double>& surfacevol () { return state_blackoil_.surfacevol(); }
std::vector<double>& facepressure() { return state_blackoil_.facepressure(); }
std::vector<double>& faceflux () { return state_blackoil_.faceflux(); }
std::vector<double>& saturation () { return state_blackoil_.saturation(); }
std::vector<double>& gasoilratio () { return state_blackoil_.gasoilratio(); }
std::vector<double>& rv () { return state_blackoil_.rv(); }
std::vector<double>& concentration() { return concentration_; }
std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& pressure () const { return state_blackoil_.pressure(); }
const std::vector<double>& temperature () const { return state_blackoil_.temperature(); }
const std::vector<double>& surfacevol () const { return state_blackoil_.surfacevol(); }
const std::vector<double>& facepressure() const { return state_blackoil_.facepressure(); }
const std::vector<double>& faceflux () const { return state_blackoil_.faceflux(); }
const std::vector<double>& saturation () const { return state_blackoil_.saturation(); }
const std::vector<double>& gasoilratio() const { return state_blackoil_.gasoilratio(); }
const std::vector<double>& rv () const { return state_blackoil_.rv(); }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& maxconcentration() const { return cmax_; }
BlackoilState& blackoilState() { return state_blackoil_; }
const BlackoilState& blackoilState() const { return state_blackoil_; }
private:
BlackoilState state_blackoil_;
std::vector<double> concentration_;
std::vector<double> cmax_;
};

52
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
View File

@ -25,7 +25,6 @@
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/fullyimplicit/PolymerBlackoilOutputWriter.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
@ -79,31 +78,29 @@
namespace Opm
{
template<class GridT>
class SimulatorFullyImplicitBlackoilPolymer;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilPolymer<GridT> >
template<>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilPolymer>
{
typedef WellStateFullyImplicitBlackoilPolymer WellState;
typedef PolymerBlackoilState ReservoirState;
typedef PolymerBlackoilOutputWriter OutputWriter;
typedef BlackoilOutputWriter OutputWriter;
typedef UnstructuredGrid Grid;
typedef BlackoilPolymerModel<Grid> Model;
};
/// Class collecting all necessary components for a two-phase simulation.
template<class GridT>
/// Class collecting all necessary components for a blackoil simulation with polymer
/// injection.
class SimulatorFullyImplicitBlackoilPolymer
: public SimulatorBase<GridT, SimulatorFullyImplicitBlackoilPolymer<GridT> >
: public SimulatorBase<SimulatorFullyImplicitBlackoilPolymer >
{
typedef SimulatorFullyImplicitBlackoilPolymer<GridT> ThisType;
typedef SimulatorBase<GridT, ThisType> BaseType;
typedef SimulatorFullyImplicitBlackoilPolymer ThisType;
typedef SimulatorBase<ThisType> BaseType;
public:
/// \brief The type of the grid that we use.
typedef GridT Grid;
SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
const Grid& grid,
const typename BaseType::Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
@ -114,21 +111,25 @@ namespace Opm
const bool vapoil,
const bool polymer,
std::shared_ptr<EclipseState> eclipse_state,
PolymerBlackoilOutputWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face);
#if 0
SimulatorReport run(SimulatorTimer& timer,
PolymerBlackoilState& state);
#endif
void computeRESV(const std::size_t step,
const Wells* wells,
const PolymerBlackoilState& x,
WellState& xw)
std::shared_ptr<Model> createModel(const typename Model::ModelParameters &modelParams,
const Wells* wells)
{
BaseType::computeRESV(step, wells, x.blackoilState(), xw);
return std::make_shared<Model>(modelParams,
BaseType::grid_,
BaseType::props_,
BaseType::geo_,
BaseType::rock_comp_props_,
polymer_props_,
wells,
BaseType::solver_,
BaseType::has_disgas_,
BaseType::has_vapoil_,
has_polymer_,
BaseType::terminal_output_);
}
void handleAdditionalWellInflow(SimulatorTimer& timer,
@ -164,4 +165,5 @@ namespace Opm
} // namespace Opm
#include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED

165
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
View File

@ -21,8 +21,7 @@
namespace Opm
{
template<class T>
SimulatorFullyImplicitBlackoilPolymer<T>::
SimulatorFullyImplicitBlackoilPolymer::
SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
@ -35,7 +34,7 @@ namespace Opm
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
PolymerBlackoilOutputWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face)
: BaseType(param,
@ -55,164 +54,4 @@ namespace Opm
, deck_(deck)
{
}
#if 0
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::run(SimulatorTimer& timer,
PolymerBlackoilState& state)
{
WellStateFullyImplicitBlackoilPolymer prev_well_state;
// Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
std::string tstep_filename = BaseType::output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
typedef T Grid;
typedef BlackoilPolymerModel<Grid> Model;
typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( BaseType::param_ );
typedef NewtonSolver<Model> Solver;
typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( BaseType::param_ );
//adaptive time stepping
// std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
// if( BaseType::param_.getDefault("timestep.adaptive", bool(false) ) )
// {
// adaptiveTimeStepping.reset( new AdaptiveTimeStepping( BaseType::param_ ) );
// }
// init output writer
BaseType::output_writer_.writeInit( timer );
std::string restorefilename = BaseType::param_.getDefault("restorefile", std::string("") );
if( ! restorefilename.empty() )
{
// -1 means that we'll take the last report step that was written
const int desiredRestoreStep = BaseType::param_.getDefault("restorestep", int(-1) );
BaseType::output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep );
}
unsigned int totalNewtonIterations = 0;
unsigned int totalLinearIterations = 0;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
step_timer.start();
if ( BaseType::terminal_output_ )
{
timer.report(std::cout);
}
// Create wells and well state.
WellsManager wells_manager(BaseType::eclipse_state_,
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(BaseType::grid_),
Opm::UgGridHelpers::globalCell(BaseType::grid_),
Opm::UgGridHelpers::cartDims(BaseType::grid_),
Opm::UgGridHelpers::dimensions(BaseType::grid_),
Opm::UgGridHelpers::cell2Faces(BaseType::grid_),
Opm::UgGridHelpers::beginFaceCentroids(BaseType::grid_),
BaseType::props_.permeability());
const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoilPolymer well_state;
well_state.init(wells, state.blackoilState(), prev_well_state);
// compute polymer inflow
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
well_state.polymerInflow() = polymer_inflow_c;
// write simulation state at the report stage
BaseType::output_writer_.writeTimeStep( timer, state.blackoilState(), well_state );
// Max oil saturation (for VPPARS), hysteresis update.
BaseType::props_.updateSatOilMax(state.saturation());
BaseType::props_.updateSatHyst(state.saturation(), BaseType::allcells_);
// Compute reservoir volumes for RESV controls.
BaseType::computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state);
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
Model model(modelParams, BaseType::grid_, BaseType::props_, BaseType::geo_, BaseType::rock_comp_props_, polymer_props_, wells, BaseType::solver_, BaseType::has_disgas_, BaseType::has_vapoil_, has_polymer_, BaseType::terminal_output_);
if (!BaseType::threshold_pressures_by_face_.empty()) {
model.setThresholdPressures(BaseType::threshold_pressures_by_face_);
}
Solver solver(solverParams, model);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge
//
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
// if( adaptiveTimeStepping ) {
// adaptiveTimeStepping->step( timer, solver, state, well_state, BaseType::output_writer_ );
// } else {
// solve for complete report step
solver.step(timer.currentStepLength(), state, well_state);
// }
// take time that was used to solve system for this reportStep
solver_timer.stop();
// accumulate the number of Newton and Linear Iterations
totalNewtonIterations += solver.newtonIterations();
totalLinearIterations += solver.linearIterations();
// Report timing.
const double st = solver_timer.secsSinceStart();
if ( BaseType::terminal_output_ )
{
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
}
stime += st;
if ( BaseType::output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
// Increment timer, remember well state.
++timer;
prev_well_state = well_state;
}
// Write final simulation state.
BaseType::output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state );
// Stop timer and create timing report
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
report.total_newton_iterations = totalNewtonIterations;
report.total_linear_iterations = totalLinearIterations;
return report;
}
#endif
} // namespace Opm

0
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

21
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp
View File

@ -30,7 +30,7 @@
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
#include <opm/polymer/fullyimplicit/PolymerBlackoilOutputWriter.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
@ -68,24 +68,25 @@
namespace Opm
{
template <class GridT>
class SimulatorFullyImplicitCompressiblePolymer;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitCompressiblePolymer<GridT> >
template <>
struct SimulatorTraits<SimulatorFullyImplicitCompressiblePolymer>
{
typedef PolymerBlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef PolymerBlackoilOutputWriter OutputWriter;
typedef BlackoilOutputWriter OutputWriter;
typedef UnstructuredGrid Grid;
#warning TODO: the 2p polymer solver does not yet adhere to The New Order!
typedef BlackoilPolymerModel<Grid> Model;
};
/// Class collecting all necessary components for a two-phase simulation.
template <class GridT>
class SimulatorFullyImplicitCompressiblePolymer
: public SimulatorBase<GridT, SimulatorFullyImplicitCompressiblePolymer<GridT> >
: public SimulatorBase<SimulatorFullyImplicitCompressiblePolymer>
{
typedef SimulatorFullyImplicitCompressiblePolymer<GridT> ThisType;
typedef SimulatorBase<GridT, ThisType> BaseType;
typedef SimulatorFullyImplicitCompressiblePolymer ThisType;
typedef SimulatorBase<ThisType> BaseType;
public:
/// Initialise from parameters and objects to observe.
@ -96,7 +97,7 @@ namespace Opm
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
PolymerBlackoilOutputWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity);

338
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp Normal file
View File

@ -0,0 +1,338 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_IMPL_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_IMPL_HEADER_INCLUDED
namespace Opm
{
namespace
{
static void outputStateVtk(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/vtk_files";
boost::filesystem::path fpath(vtkfilename.str());
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
}
static void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
dm["surfvolume"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
#if 0
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
std::ofstream os(fname.c_str());
if (!os) {
OPM_THROW(std::runtime_error, "Failed to open " << fname);
}
watercut.write(os);
}
#endif
}
/// Class collecting all necessary components for a two-phase simulation.
SimulatorFullyImplicitCompressiblePolymer::
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
: BaseType(param,
grid,
geo,
props,
rock_comp_props,
linsolver,
gravity,
/*disgas=*/false,
/*vapoil=*/false,
eclipse_state,
output_writer,
/*threshold_pressures_by_face=*/std::vector<double>())
, deck_(deck)
, polymer_props_(polymer_props)
{
// For output.
output_ = param.getDefault("output", true);
if (output_) {
output_vtk_ = param.getDefault("output_vtk", true);
output_dir_ = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir_);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
output_interval_ = param.getDefault("output_interval", 1);
}
// Well control related init.
check_well_controls_ = param.getDefault("check_well_controls", false);
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Misc init.
const int num_cells = grid.number_of_cells;
BaseType::allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
BaseType::allcells_[cell] = cell;
}
}
SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer,
typename BaseType::ReservoirState& state)
{
WellStateFullyImplicitBlackoil prev_well_state;
// Initialisation.
std::vector<double> porevol;
if (BaseType::rock_comp_props_ && BaseType::rock_comp_props_->isActive()) {
computePorevolume(BaseType::grid_,
BaseType::props_.porosity(),
*BaseType::rock_comp_props_,
state.pressure(),
porevol);
} else {
computePorevolume(BaseType::grid_,
BaseType::props_.porosity(),
porevol);
}
std::vector<double> initial_porevol = porevol;
std::vector<double> polymer_inflow_c(BaseType::grid_.number_of_cells);
// Main simulation loop.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
std::string tstep_filename = output_dir_ + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
//Main simulation loop.
while (!timer.done()) {
#if 0
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
if (wells_) {
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
}
std::fstream tstep_os;
if (output_) {
std::string filename = output_dir_ + "/step_timing.param";
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
}
#endif
// Report timestep and (optionally) write state to disk.
step_timer.start();
timer.report(std::cout);
WellsManager wells_manager(BaseType::eclipse_state_,
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(BaseType::grid_),
Opm::UgGridHelpers::globalCell(BaseType::grid_),
Opm::UgGridHelpers::cartDims(BaseType::grid_),
Opm::UgGridHelpers::dimensions(BaseType::grid_),
Opm::UgGridHelpers::cell2Faces(BaseType::grid_),
Opm::UgGridHelpers::beginFaceCentroids(BaseType::grid_),
BaseType::props_.permeability());
const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoil well_state;
well_state.init(wells, state, prev_well_state);
//Compute polymer inflow.
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
if (output_vtk_) {
outputStateVtk(BaseType::grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(BaseType::grid_, state, timer.currentStepNum(), output_dir_);
}
if (output_) {
if (timer.currentStepNum() == 0) {
BaseType::output_writer_.writeInit(timer);
}
BaseType::output_writer_.writeTimeStep(timer, state, well_state);
}
// Run solver.
solver_timer.start();
FullyImplicitCompressiblePolymerSolver solver(BaseType::grid_, BaseType::props_, BaseType::geo_, BaseType::rock_comp_props_, polymer_props_, *wells_manager.c_wells(), BaseType::solver_);
solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
// Stop timer and report.
solver_timer.stop();
const double st = solver_timer.secsSinceStart();
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
stime += st;
// Update pore volumes if rock is compressible.
if (BaseType::rock_comp_props_ && BaseType::rock_comp_props_->isActive()) {
initial_porevol = porevol;
computePorevolume(BaseType::grid_, BaseType::props_.porosity(), *BaseType::rock_comp_props_, state.pressure(), porevol);
}
/*
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double polyinj = 0;
double polyprod = 0;
Opm::computeInjectedProduced(props_, polymer_props_,
state,
transport_src, polymer_inflow_c, timer.currentStepLength(),
injected, produced,
polyinj, polyprod);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
std::cout.precision(5);
const int width = 18;
std::cout << "\nMass balance report.\n";
std::cout << " Injected reservoir volumes: "
<< std::setw(width) << injected[0]
<< std::setw(width) << injected[1] << std::endl;
std::cout << " Produced reservoir volumes: "
<< std::setw(width) << produced[0]
<< std::setw(width) << produced[1] << std::endl;
std::cout << " Total inj reservoir volumes: "
<< std::setw(width) << tot_injected[0]
<< std::setw(width) << tot_injected[1] << std::endl;
std::cout << " Total prod reservoir volumes: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1] << std::endl;
*/
if (output_) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
++timer;
prev_well_state = well_state;
}
// Write final simulation state.
if (output_) {
if (output_vtk_) {
outputStateVtk(BaseType::grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(BaseType::grid_, state, timer.currentStepNum(), output_dir_);
BaseType::output_writer_.writeTimeStep(timer, state, prev_well_state);
}
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
return report;
}
} // namespace Opm
#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED