Merge pull request #1081 from atgeirr/restart-no-convert-call

Improve restarting
This commit is contained in:
Atgeirr Flø Rasmussen 2017-03-03 14:09:58 +01:00 committed by GitHub
commit 4bb064b329
8 changed files with 146 additions and 33 deletions

View File

@ -75,10 +75,18 @@ std::vector< double >& stripe( const std::vector< double >& v,
data::Solution simToSolution( const SimulationDataContainer& reservoir,
const bool use_si_units,
PhaseUsage phases ) {
// Set up unit system to use to suppress conversion if use_si_units is true.
const UnitSystem::measure press_unit = use_si_units ? UnitSystem::measure::identity : UnitSystem::measure::pressure;
const UnitSystem::measure temp_unit = use_si_units ? UnitSystem::measure::identity : UnitSystem::measure::temperature;
const UnitSystem::measure rs_unit = use_si_units ? UnitSystem::measure::identity : UnitSystem::measure::gas_oil_ratio;
const UnitSystem::measure rv_unit = use_si_units ? UnitSystem::measure::identity : UnitSystem::measure::oil_gas_ratio;
data::Solution sol;
sol.insert( "PRESSURE", UnitSystem::measure::pressure, reservoir.pressure() , data::TargetType::RESTART_SOLUTION);
sol.insert( "TEMP" , UnitSystem::measure::temperature, reservoir.temperature() , data::TargetType::RESTART_SOLUTION );
sol.insert( "PRESSURE", press_unit, reservoir.pressure() , data::TargetType::RESTART_SOLUTION);
sol.insert( "TEMP" , temp_unit, reservoir.temperature() , data::TargetType::RESTART_SOLUTION );
const auto ph = reservoir.numPhases();
const auto& sat = reservoir.saturation();
@ -95,11 +103,11 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir,
}
if( reservoir.hasCellData( BlackoilState::GASOILRATIO ) ) {
sol.insert( "RS", UnitSystem::measure::gas_oil_ratio, reservoir.getCellData( BlackoilState::GASOILRATIO ) , data::TargetType::RESTART_SOLUTION );
sol.insert( "RS", rs_unit, reservoir.getCellData( BlackoilState::GASOILRATIO ) , data::TargetType::RESTART_SOLUTION );
}
if( reservoir.hasCellData( BlackoilState::RV ) ) {
sol.insert( "RV", UnitSystem::measure::oil_gas_ratio, reservoir.getCellData( BlackoilState::RV ) , data::TargetType::RESTART_SOLUTION );
sol.insert( "RV", rv_unit, reservoir.getCellData( BlackoilState::RV ) , data::TargetType::RESTART_SOLUTION );
}
if (reservoir.hasCellData( BlackoilSolventState::SSOL)) {

View File

@ -51,7 +51,10 @@ namespace Opm {
/// Returns Solution with the following fields:
/// PRESSURE, TEMP (unconditionally)
/// SWAT, SGAS, RS, RV, SSOL (if appropriate fields present in input)
/// If use_si_units is true, the fields will have the measure UnitSystem::measure::identity,
/// and therefore *not* be converted to customary units (depending on family) upon output.
data::Solution simToSolution( const SimulationDataContainer& reservoir,
const bool use_si_units,
PhaseUsage phases );
/// Copies the following fields from sol into state (all conditionally):

View File

@ -89,10 +89,10 @@ namespace Opm
{
WellState prev_well_state;
ExtraData extra;
if (output_writer_.isRestart()) {
// This is a restart, populate WellState and ReservoirState state objects from restart file
output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state);
output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state, extra);
initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid_), has_disgas_, has_vapoil_);
}
@ -117,6 +117,9 @@ namespace Opm
} else {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
}
if (output_writer_.isRestart()) {
adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
}
}
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
@ -188,7 +191,7 @@ namespace Opm
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state, {} );
report.output_write_time += perfTimer.stop();
}

View File

@ -144,9 +144,10 @@ public:
{
WellState prev_well_state;
ExtraData extra;
if (output_writer_.isRestart()) {
// This is a restart, populate WellState and ReservoirState state objects from restart file
output_writer_.initFromRestartFile(props_.phaseUsage(), grid(), state, prev_well_state);
output_writer_.initFromRestartFile(props_.phaseUsage(), grid(), state, prev_well_state, extra);
initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
}
@ -170,6 +171,10 @@ public:
} else {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
}
if (output_writer_.isRestart()) {
adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
}
}
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
@ -354,7 +359,8 @@ public:
// write simulation state at the report stage
Dune::Timer perfTimer;
perfTimer.start();
output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
output_writer_.writeTimeStep( timer, state, well_state, solver->model(), false, nextstep );
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;

View File

@ -88,6 +88,19 @@ namespace Opm
unsigned int totalLinearIterations = 0;
DynamicListEconLimited dynamic_list_econ_limited;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
//Get compressed cell fipnum.
std::vector<int> fipnum(AutoDiffGrid::numCells(grid_));
if (fipnum_global.empty()) {
std::fill(fipnum.begin(), fipnum.end(), 0);
} else {
for (size_t c = 0; c < fipnum.size(); ++c) {
fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
}
}
std::vector<std::vector<double> > OOIP;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@ -132,7 +145,7 @@ namespace Opm
if (timer.initialStep()) {
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state, {} );
}
// Max oil saturation (for VPPARS), hysteresis update.
@ -147,6 +160,13 @@ namespace Opm
auto solver = createSolver(well_model);
// Compute orignal FIP;
if (!ooip_computed) {
OOIP = solver->computeFluidInPlace(state, fipnum);
Base::FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
ooip_computed = true;
}
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are too large for the solver to converge
//
@ -170,6 +190,25 @@ namespace Opm
// Report timing.
const double st = solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<std::vector<double> > COIP;
COIP = solver->computeFluidInPlace(state, fipnum);
std::vector<double> OOIP_totals = Base::FIPTotals(OOIP, state);
std::vector<double> COIP_totals = Base::FIPTotals(COIP, state);
//Convert to correct units
Base::FIPUnitConvert(eclipse_state_->getUnits(), COIP);
Base::FIPUnitConvert(eclipse_state_->getUnits(), OOIP_totals);
Base::FIPUnitConvert(eclipse_state_->getUnits(), COIP_totals);
if ( terminal_output_ )
{
Base::outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
Base::outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
}
}
if ( terminal_output_ )
{
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;

View File

@ -203,6 +203,7 @@ namespace Opm
const SimulationDataContainer state_;
const WellStateFullyImplicitBlackoil wellState_;
data::Solution simProps_;
std::map<std::string, std::vector<double>> extraData_;
const bool substep_;
explicit WriterCall( BlackoilOutputWriter& writer,
@ -210,12 +211,14 @@ namespace Opm
const SimulationDataContainer& state,
const WellStateFullyImplicitBlackoil& wellState,
const data::Solution& simProps,
const std::map<std::string, std::vector<double>>& extraData,
bool substep )
: writer_( writer ),
timer_( timer.clone() ),
state_( state ),
wellState_( wellState ),
simProps_( simProps ),
extraData_( extraData ),
substep_( substep )
{
}
@ -224,7 +227,7 @@ namespace Opm
void run ()
{
// write data
writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ );
writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, extraData_, substep_ );
}
};
}
@ -238,15 +241,16 @@ namespace Opm
const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellStateFullyImplicitBlackoil& localWellState,
const std::map<std::string, std::vector<double>>& extraData,
bool substep)
{
data::Solution localCellData{};
if( output_ )
{
localCellData = simToSolution(localState, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
localCellData = simToSolution(localState, restart_double_si_, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
}
writeTimeStepWithCellProperties(timer, localState, localCellData ,
localWellState, substep);
localWellState, extraData, substep);
}
@ -260,6 +264,7 @@ namespace Opm
const SimulationDataContainer& localState,
const data::Solution& localCellData,
const WellStateFullyImplicitBlackoil& localWellState,
const std::map<std::string, std::vector<double>>& extraData,
bool substep)
{
// VTK output (is parallel if grid is parallel)
@ -281,6 +286,7 @@ namespace Opm
isIORank = parallelOutput_->collectToIORank( localState, localWellState,
localCellData,
wellStateStepNumber );
// Note that at this point the extraData are assumed to be global, i.e. identical across all processes.
}
const data::Solution& cellData = ( parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalCellData() : localCellData;
@ -294,12 +300,12 @@ namespace Opm
{
if( asyncOutput_ ) {
// dispatch the write call to the extra thread
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, cellData, substep ) );
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, cellData, extraData, substep ) );
}
else {
// just write the data to disk
try {
writeTimeStepSerial( timer, state, wellState, cellData, substep );
writeTimeStepSerial( timer, state, wellState, cellData, extraData, substep );
} catch (std::runtime_error& msg) {
err = 1;
emsg = msg.what();
@ -329,6 +335,7 @@ namespace Opm
const SimulationDataContainer& state,
const WellStateFullyImplicitBlackoil& wellState,
const data::Solution& simProps,
const std::map<std::string, std::vector<double>>& extraData,
bool substep)
{
// Matlab output
@ -348,7 +355,9 @@ namespace Opm
substep,
timer.simulationTimeElapsed(),
simProps,
wellState.report(phaseUsage_));
wellState.report(phaseUsage_),
extraData,
restart_double_si_);
}
}
@ -416,7 +425,7 @@ namespace Opm
// No per cell data is written for restore steps, but will be
// for subsequent steps, when we have started simulating
writeTimeStepWithoutCellProperties( timer, state, wellState );
writeTimeStepWithoutCellProperties( timer, state, wellState, {} );
// some output
std::cout << "Restored step " << timer.reportStepNum() << " at day "

View File

@ -50,6 +50,8 @@
#include <iomanip>
#include <fstream>
#include <thread>
#include <map>
#include <set>
#include <boost/filesystem.hpp>
@ -197,6 +199,14 @@ namespace Opm
const Grid& grid_;
};
/// Extra data to read/write for OPM restarting
struct ExtraData
{
double suggested_step = -1.0;
};
/** \brief Wrapper class for VTK, Matlab, and ECL output. */
class BlackoilOutputWriter
{
@ -224,19 +234,21 @@ namespace Opm
const SimulationDataContainer& reservoirState,
const Opm::WellStateFullyImplicitBlackoil& wellState,
const Model& physicalModel,
bool substep = false);
bool substep = false,
const double nextstep = -1.0);
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight. This function will write all
* CellData in simProps to the file as well.
* CellData in simProps to the file as well as the extraData.
*/
void writeTimeStepWithCellProperties(
const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const data::Solution& cellData,
const Opm::WellStateFullyImplicitBlackoil& wellState,
const std::map<std::string, std::vector<double>>& extraData,
bool substep = false);
/*!
@ -248,6 +260,7 @@ namespace Opm
const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellStateFullyImplicitBlackoil& wellState,
const std::map<std::string, std::vector<double>>& extraData,
bool substep = false);
/*!
@ -259,6 +272,7 @@ namespace Opm
const SimulationDataContainer& reservoirState,
const Opm::WellStateFullyImplicitBlackoil& wellState,
const data::Solution& simProps,
const std::map<std::string, std::vector<double>>& extraData,
bool substep);
/** \brief return output directory */
@ -280,11 +294,12 @@ namespace Opm
const int desiredReportStep);
template <class Grid, class WellStateFullyImplicitBlackOel>
template <class Grid, class WellState>
void initFromRestartFile(const PhaseUsage& phaseUsage,
const Grid& grid,
SimulationDataContainer& simulatorstate,
WellStateFullyImplicitBlackOel& wellstate);
WellState& wellstate,
ExtraData& extra);
bool isRestart() const;
@ -297,6 +312,7 @@ namespace Opm
// Parameters for output.
const std::string outputDir_;
const int output_interval_;
const bool restart_double_si_;
int lastBackupReportStep_;
@ -328,6 +344,7 @@ namespace Opm
parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, phaseUsage ) : 0 ),
outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ),
output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ),
restart_double_si_( output_ ? param.getDefault("restart_double_si", false) : false ),
lastBackupReportStep_( -1 ),
phaseUsage_( phaseUsage ),
eclipseState_(eclipseState),
@ -401,7 +418,8 @@ namespace Opm
initFromRestartFile( const PhaseUsage& phaseUsage,
const Grid& grid,
SimulationDataContainer& simulatorstate,
WellState& wellstate)
WellState& wellstate,
ExtraData& extra )
{
std::map<std::string, UnitSystem::measure> solution_keys {{"PRESSURE" , UnitSystem::measure::pressure},
{"SWAT" , UnitSystem::measure::identity},
@ -409,6 +427,14 @@ namespace Opm
{"TEMP" , UnitSystem::measure::temperature},
{"RS" , UnitSystem::measure::gas_oil_ratio},
{"RV" , UnitSystem::measure::oil_gas_ratio}};
std::set<std::string> extra_keys {"OPMEXTRA"};
if (restart_double_si_) {
// Avoid any unit conversions, treat restart input as SI units.
for (auto& elem : solution_keys) {
elem.second = UnitSystem::measure::identity;
}
}
// gives a dummy dynamic_list_econ_limited
DynamicListEconLimited dummy_list_econ_limited;
@ -431,10 +457,19 @@ namespace Opm
const Wells* wells = wellsmanager.c_wells();
wellstate.resize(wells, simulatorstate, phaseUsage ); //Resize for restart step
auto state = eclIO_->loadRestart(solution_keys);
auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
solutionToSim( state.first, phaseUsage, simulatorstate );
wellsToState( state.second, phaseUsage, wellstate );
solutionToSim( restart_values.solution, phaseUsage, simulatorstate );
wellsToState( restart_values.wells, phaseUsage, wellstate );
const auto opmextra_iter = restart_values.extra.find("OPMEXTRA");
if (opmextra_iter != restart_values.extra.end()) {
std::vector<double> opmextra = opmextra_iter->second;
assert(opmextra.size() == 1);
extra.suggested_step = opmextra[0];
} else {
OPM_THROW(std::runtime_error, "Cannot restart, restart data is missing OPMEXTRA field.");
}
}
@ -899,13 +934,15 @@ namespace Opm
const SimulationDataContainer& localState,
const WellStateFullyImplicitBlackoil& localWellState,
const Model& physicalModel,
bool substep)
bool substep,
const double nextstep)
{
data::Solution localCellData{};
const RestartConfig& restartConfig = eclipseState_.getRestartConfig();
const SummaryConfig& summaryConfig = eclipseState_.getSummaryConfig();
const int reportStepNum = timer.reportStepNum();
bool logMessages = output_ && parallelOutput_->isIORank();
std::map<std::string, std::vector<double>> extraData;
if( output_ )
{
@ -917,16 +954,20 @@ namespace Opm
SimulationDataContainer sd =
detail::convertToSimulationDataContainer( physicalModel.getSimulatorData(), localState, phaseUsage_ );
localCellData = simToSolution( sd, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
localCellData = simToSolution( sd, restart_double_si_, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
detail::getRestartData( localCellData, std::move(sd), phaseUsage_, physicalModel,
restartConfig, reportStepNum, logMessages );
// sd will be invalid after getRestartData has been called
}
detail::getSummaryData( localCellData, phaseUsage_, physicalModel, summaryConfig );
assert(!localCellData.empty());
// Add suggested next timestep to extra data.
extraData["OPMEXTRA"] = std::vector<double>(1, nextstep);
}
writeTimeStepWithCellProperties(timer, localState, localCellData, localWellState, substep);
writeTimeStepWithCellProperties(timer, localState, localCellData, localWellState, extraData, substep);
}
}
#endif

View File

@ -84,6 +84,10 @@ namespace Opm {
Output& outputWriter,
const std::vector<int>* fipnum = nullptr);
double suggestedNextStep() const { return suggested_next_timestep_; }
void setSuggestedNextStep(const double x) { suggested_next_timestep_ = x; }
protected:
template <class Solver, class State, class WellState, class Output>
SimulatorReport stepImpl( const SimulatorTimer& timer,