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

Merge pull request #179 from atgeirr/restructure-sim

Browse Source 
Restructure fully implicit blackoil simulators
This commit is contained in:
Bård Skaflestad 2014-08-13 10:23:41 +02:00
commit 6fccb4b446
4 changed files with 121 additions and 217 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
examples/sim_fibo_ad.cpp
View File

@ -163,9 +163,9 @@ try
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::ofstream outStream;
std::string output_dir;
if (output) {
// Create output directory if needed.
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
@ -175,20 +175,11 @@ try
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
std::string filename = output_dir + "/timing.param";
outStream.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
// open file to clean it. The file is appended to in SimulatorTwophase
filename = output_dir + "/step_timing.param";
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
step_os.close();
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
WellStateFullyImplicitBlackoil well_state;
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
@ -196,59 +187,28 @@ try
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of a report step.
std::cout << "\n"
<< "---------------------------------------------------------------\n"
<< "-------------- Starting report step " << reportStepIdx << " --------------\n"
<< "---------------------------------------------------------------\n"
<< "\n";
SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
*grid->c_grid(),
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
eclipseState,
outputWriter);
WellsManager wells(eclipseState,
reportStepIdx,
*grid->c_grid(),
props->permeability());
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
if (reportStepIdx == 0) {
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.
well_state.init(wells.c_wells(), state);
}
simtimer.setCurrentStepNum(reportStepIdx);
if (reportStepIdx == 0) {
outputWriter.writeInit(simtimer);
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
}
// Create and run simulator.
SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
*grid->c_grid(),
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL") );
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
++simtimer;
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
fullReport += episodeReport;
}
SimulatorReport fullReport = simulator.run(simtimer, state);
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.report(std::cout);
if (output) {
std::string filename = output_dir + "/walltime.param";
std::string filename = output_dir + "/walltime.txt";
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
fullReport.reportParam(tot_os);
warnIfUnusedParams(param);

83
examples/sim_fibo_ad_cp.cpp
View File

@ -48,6 +48,7 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
@ -130,8 +131,9 @@ try
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr newParser(new Opm::Parser() );
Opm::DeckConstPtr deck = newParser->parseFile( deck_filename );
Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck));
bool strict_parsing = param.getDefault("strict_parsing", true);
Opm::DeckConstPtr deck = newParser->parseFile(deck_filename, strict_parsing);
std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck));
// Grid init
grid.reset(new Dune::CpGrid());
@ -184,6 +186,8 @@ try
/ state.surfacevol()[c*np + pu.phase_pos[Oil]];
}
}
} else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
OPM_THROW(std::logic_error, "sim_fibo_ad_cp does not support EQUIL initialization.");
} else {
initBlackoilStateFromDeck(grid->numCells(), &(grid->globalCell())[0],
grid->numFaces(), AutoDiffGrid::faceCells(*grid),
@ -205,9 +209,9 @@ try
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::ofstream outStream;
std::string output_dir;
if (output) {
// Create output directory if needed.
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
@ -217,19 +221,13 @@ try
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
std::string filename = output_dir + "/timing.param";
outStream.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
// open file to clean it. The file is appended to in SimulatorTwophase
filename = output_dir + "/step_timing.param";
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
step_os.close();
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
WellStateFullyImplicitBlackoil well_state;
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer;
@ -238,63 +236,28 @@ try
Opm::DerivedGeology geology(*grid, *new_props, eclipseState, grav);
SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of a report step.
std::cout << "\n"
<< "---------------------------------------------------------------\n"
<< "-------------- Starting report step " << reportStepIdx << " --------------\n"
<< "---------------------------------------------------------------\n"
<< "\n";
SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
*grid,
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
eclipseState,
outputWriter);
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, Opm::UgGridHelpers::numCells(*grid),
Opm::UgGridHelpers::globalCell(*grid),
Opm::UgGridHelpers::cartDims(*grid),
Opm::UgGridHelpers::dimensions(*grid),
Opm::UgGridHelpers::beginCellCentroids(*grid),
Opm::UgGridHelpers::cell2Faces(*grid),
Opm::UgGridHelpers::beginFaceCentroids(*grid),
props->permeability());
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
if (reportStepIdx == 0) {
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.
well_state.init(wells.c_wells(), state);
}
simtimer.setCurrentStepNum(reportStepIdx);
if (reportStepIdx == 0) {
outputWriter.writeInit(simtimer);
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
}
SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
*grid,
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL") );
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
++simtimer;
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
fullReport += episodeReport;
}
SimulatorReport fullReport = simulator.run(simtimer, state);
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.report(std::cout);
if (output) {
std::string filename = output_dir + "/walltime.param";
std::string filename = output_dir + "/walltime.txt";
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
fullReport.reportParam(tot_os);
warnIfUnusedParams(param);

12
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -38,6 +38,8 @@ namespace Opm
class SimulatorTimer;
class BlackoilState;
class WellStateFullyImplicitBlackoil;
class EclipseState;
class EclipseWriter;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation.
@ -69,16 +71,19 @@ namespace Opm
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] eclipse_state
/// \param[in] output_writer
SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil );
const bool vapoil,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will
@ -88,8 +93,7 @@ namespace Opm
/// \param[in,out] well_state state of wells: bhp, perforation rates
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
BlackoilState& state);
private:
class Impl;

167
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -30,6 +30,7 @@
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
@ -65,19 +66,19 @@ namespace Opm
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
bool has_disgas,
bool has_vapoil );
bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer);
SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
BlackoilState& state);
private:
// Data.
const parameter::ParameterGroup param_;
// Parameters for output.
bool output_;
bool output_vtk_;
@ -90,14 +91,18 @@ namespace Opm
const Grid& grid_;
BlackoilPropsAdInterface& props_;
const RockCompressibility* rock_comp_props_;
WellsManager& wells_manager_;
const Wells* wells_;
const double* gravity_;
// Solvers
const DerivedGeology &geo_;
FullyImplicitBlackoilSolver<Grid> solver_;
const DerivedGeology& geo_;
NewtonIterationBlackoilInterface& solver_;
// Misc. data
std::vector<int> allcells_;
const bool has_disgas_;
const bool has_vapoil_;
// eclipse_state
std::shared_ptr<EclipseState> eclipse_state_;
// output_writer
EclipseWriter& output_writer_;
};
@ -109,14 +114,16 @@ namespace Opm
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil )
const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer)
{
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, wells_manager, linsolver, gravity, has_disgas, has_vapoil));
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil,
eclipse_state, output_writer));
}
@ -125,10 +132,9 @@ namespace Opm
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state)
BlackoilState& state)
{
return pimpl_->run(timer, state, well_state);
return pimpl_->run(timer, state);
}
@ -197,23 +203,23 @@ namespace Opm
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil)
: grid_(grid),
const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer)
: param_(param),
grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
gravity_(gravity),
geo_(geo),
solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver, has_disgas, has_vapoil)
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, */
solver_(linsolver),
has_disgas_(has_disgas),
has_vapoil_(has_vapoil),
eclipse_state_(eclipse_state),
output_writer_(output_writer)
{
// For output.
output_ = param.getDefault("output", true);
@ -245,39 +251,45 @@ namespace Opm
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state)
BlackoilState& state)
{
// Initialisation.
std::vector<double> porevol;
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
} else {
computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), porevol);
}
// const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
WellStateFullyImplicitBlackoil well_state;
// Main simulation loop.
// 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::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);
}
std::string tstep_filename = output_dir_ + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
// Main simulation loop.
while (!timer.done()) {
// Report timestep and (optionally) write state to disk.
step_timer.start();
timer.report(std::cout);
WellsManager wells_manager(eclipse_state_,
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(grid_),
Opm::UgGridHelpers::globalCell(grid_),
Opm::UgGridHelpers::cartDims(grid_),
Opm::UgGridHelpers::dimensions(grid_),
Opm::UgGridHelpers::beginCellCentroids(grid_),
Opm::UgGridHelpers::cell2Faces(grid_),
Opm::UgGridHelpers::beginFaceCentroids(grid_),
props_.permeability());
const Wells *wells = wells_manager.c_wells();
if (timer.currentStepNum() == 0) {
well_state.init(wells, state);
output_writer_.writeInit(timer);
} else {
// TODO: add a function to update the well_state here.
}
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
@ -289,58 +301,24 @@ namespace Opm
SimulatorReport sreport;
FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
// Max oil saturation
props_.updateSatOilMax(state.saturation());
// Hysteresis
props_.updateSatHyst(state.saturation(), allcells_);
// Solve pressure equation.
// if (check_well_controls_) {
// computeFractionalFlow(props_, allcells_,
// state.pressure(), state.surfacevol(), state.saturation(),
// fractional_flows);
// wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
// }
// Run solver.
solver_timer.start();
solver.step(timer.currentStepLength(), state, well_state);
bool well_control_passed = !check_well_controls_;
int well_control_iteration = 0;
do {
// Run solver.
solver_timer.start();
std::vector<double> initial_pressure = state.pressure();
solver_.step(timer.currentStepLength(), state, well_state);
// Stop timer and report.
solver_timer.stop();
const double st = solver_timer.secsSinceStart();
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
// 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;
sreport.pressure_time = st;
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_state.wellRates(), well_state.wellRates());
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
if (rock_comp_props_ && rock_comp_props_->isActive()) {
initial_porevol = porevol;
computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
}
stime += st;
sreport.pressure_time = st;
sreport.total_time = step_timer.secsSinceStart();
if (output_) {
@ -351,12 +329,11 @@ namespace Opm
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
tstep_os.close();
}
// advance to next timestep before reporting at this location
// ++timer; // Commented out since this has temporarily moved to the main() function.
break; // this is a temporary measure
output_writer_.writeTimeStep(timer, state, well_state.basicWellState());
++timer;
}
total_timer.stop();