Initial version for outputting cell data

This commit is contained in:
babrodtk 2016-09-01 11:38:04 +02:00
parent 4d2d004c43
commit 35bed24465
7 changed files with 239 additions and 81 deletions

View File

@ -268,6 +268,17 @@ namespace Opm {
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (pu.phase_used[phase]) {
const int pos = pu.phase_pos[phase];
return rq_[pos].b;
}
else {
return ADB::null();
}
}
protected:
// --------- Types and enums ---------

View File

@ -260,6 +260,10 @@ namespace Opm {
}
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
return transport_model_->getReciprocalFormationVolumeFactor(phase);
}
protected:
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;

View File

@ -176,7 +176,10 @@ namespace Opm
// write the inital state at the report stage
if (timer.initialStep()) {
output_writer_.writeTimeStep( timer, state, well_state );
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
std::vector<data::CellData> noData;
output_writer_.writeTimeStep( timer, state, well_state, noData );
}
// Max oil saturation (for VPPARS), hysteresis update.
@ -191,7 +194,7 @@ namespace Opm
const WellModel well_model(wells);
auto solver = asImpl().createSolver(well_model);
std::unique_ptr<Solver> solver = asImpl().createSolver(well_model);
// Compute orignal FIP;
if (!ooip_computed) {
@ -296,7 +299,8 @@ namespace Opm
++timer;
// write simulation state at the report stage
output_writer_.writeTimeStep( timer, state, well_state );
const auto& physicalModel = solver->model();
output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
prev_well_state = well_state;
// The well potentials are only computed if they are needed

View File

@ -125,7 +125,8 @@ namespace Opm
// write the inital state at the report stage
if (timer.initialStep()) {
output_writer_.writeTimeStep( timer, state, well_state );
std::vector<data::CellData> noData;
output_writer_.writeTimeStep( timer, state, well_state, noData );
}
// Max oil saturation (for VPPARS), hysteresis update.
@ -179,8 +180,10 @@ namespace Opm
// Increment timer, remember well state.
++timer;
// write simulation state at the report stage
output_writer_.writeTimeStep( timer, state, well_state );
const auto& physicalModel = solver->model();
output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
prev_well_state = well_state;
}

View File

@ -42,6 +42,11 @@
#include <boost/filesystem.hpp>
//For OutputWriterHelper
#include <map>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#ifdef HAVE_OPM_GRID
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
@ -52,6 +57,17 @@ namespace Opm
{
void outputStateVtk(const UnstructuredGrid& grid,
const SimulationDataContainer& state,
const int step,
@ -247,87 +263,14 @@ namespace Opm
namespace detail {
struct WriterCall : public ThreadHandle :: ObjectInterface
{
BlackoilOutputWriter& writer_;
std::unique_ptr< SimulatorTimerInterface > timer_;
const SimulationDataContainer state_;
const WellState wellState_;
const bool substep_;
explicit WriterCall( BlackoilOutputWriter& writer,
const SimulatorTimerInterface& timer,
const SimulationDataContainer& state,
const WellState& wellState,
bool substep )
: writer_( writer ),
timer_( timer.clone() ),
state_( state ),
wellState_( wellState ),
substep_( substep )
{
}
// callback to writer's serial writeTimeStep method
void run ()
{
// write data
writer_.writeTimeStepSerial( *timer_, state_, wellState_, substep_ );
}
};
}
void
BlackoilOutputWriter::
writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellState& localWellState,
bool substep)
{
// VTK output (is parallel if grid is parallel)
if( vtkWriter_ ) {
vtkWriter_->writeTimeStep( timer, localState, localWellState, false );
}
bool isIORank = output_ ;
if( parallelOutput_ && parallelOutput_->isParallel() )
{
// If this is not the initial write and no substep, then the well
// state used in the computation is actually the one of the last
// step. We need that well state for the gathering. Otherwise
// It an exception with a message like "global state does not
// contain well ..." might be thrown.
int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
(timer.reportStepNum() - 1) : timer.reportStepNum();
// collect all solutions to I/O rank
isIORank = parallelOutput_->collectToIORank( localState, localWellState, wellStateStepNumber );
}
const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
// serial output is only done on I/O rank
if( isIORank )
{
if( asyncOutput_ ) {
// dispatch the write call to the extra thread
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, substep ) );
}
else {
// just write the data to disk
writeTimeStepSerial( timer, state, wellState, substep );
}
}
}
void
BlackoilOutputWriter::
writeTimeStepSerial(const SimulatorTimerInterface& timer,
const SimulationDataContainer& state,
const WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep)
{
// Matlab output
@ -342,7 +285,6 @@ namespace Opm
if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
} else {
std::vector<data::CellData> simProps;
/*
The simProps vector can be passed to the writeTimestep routine
to add more properties to the restart file. Examples of the
@ -439,7 +381,11 @@ namespace Opm
restorefile >> state;
restorefile >> wellState;
writeTimeStep( timer, state, wellState );
// FIXME: We this should optimally have the proper per cell data to dump
// Right now it will not dump any per cell data until we start simulating
std::vector<data::CellData> noData;
writeTimeStep( timer, state, wellState, noData );
// some output
std::cout << "Restored step " << timer.reportStepNum() << " at day "
<< unit::convert::to(timer.simulationTimeElapsed(),unit::day) << std::endl;

View File

@ -215,15 +215,18 @@ namespace Opm
const double* permeability );
/** \copydoc Opm::OutputWriter::writeTimeStep */
template<class Model>
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
const Model& physicalModel,
bool substep = false);
/** \copydoc Opm::OutputWriter::writeTimeStep */
void writeTimeStepSerial(const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep);
/** \brief return output directory */
@ -371,5 +374,181 @@ namespace Opm
}
namespace detail {
struct WriterCall : public ThreadHandle :: ObjectInterface
{
BlackoilOutputWriter& writer_;
std::unique_ptr< SimulatorTimerInterface > timer_;
const SimulationDataContainer state_;
const WellState wellState_;
std::vector<data::CellData> simProps_;
const bool substep_;
explicit WriterCall( BlackoilOutputWriter& writer,
const SimulatorTimerInterface& timer,
const SimulationDataContainer& state,
const WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep )
: writer_( writer ),
timer_( timer.clone() ),
state_( state ),
wellState_( wellState ),
simProps_( simProps ),
substep_( substep )
{
}
// callback to writer's serial writeTimeStep method
void run ()
{
// write data
writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ );
}
};
template<class Model>
std::vector<data::CellData> getCellData(const Model& model,
const RestartConfig& restartConfig,
const int reportStepNum) {
std::vector<data::CellData> simProps;
std::map<const char*, int> outKeywords {
{"ALLPROPS", 0},
{"BG", 0},
{"BO", 0},
{"BW", 0},
{"CONV", 0},
{"DEN", 0},
{"KRG", 0},
{"KRO", 0},
{"KRW", 0},
{"RVSAT", 0},
{"RSSAT", 0},
{"NORST", 0},
{"PBPD", 0},
{"VISC", 0}
};
//Get the value of each of the keys
for (auto& keyValue : outKeywords) {
keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
}
//Postprocess some of the special keys
if (outKeywords["ALLPROPS"] > 0) {
//ALLPROPS implies KRO,KRW,KRG,xxx_DEN,xxx_VISC,BG,BO (xxx= OIL,GAS,WAT)
outKeywords["BG"] = std::max(outKeywords["BG"], 1);
outKeywords["BO"] = std::max(outKeywords["BO"], 1);
outKeywords["BW"] = std::max(outKeywords["BW"], 1);
outKeywords["KRG"] = std::max(outKeywords["KRG"], 1);
outKeywords["KRO"] = std::max(outKeywords["KRO"], 1);
outKeywords["KRW"] = std::max(outKeywords["KRW"], 1);
outKeywords["DEN"] = std::max(outKeywords["DEN"], 1);
outKeywords["VISC"] = std::max(outKeywords["VISC"], 1);
}
if (outKeywords["BW"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Aqua);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBW", Opm::UnitSystem::measure::volume, b);
}
if (outKeywords["BO"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Liquid);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBO", Opm::UnitSystem::measure::volume, b);
}
if (outKeywords["BG"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Vapour);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBG", Opm::UnitSystem::measure::volume, b);
}
return simProps;
}
/**
* Template specialization to print raw cell data
*/
template<>
inline
std::vector<data::CellData> getCellData<std::vector<data::CellData> >(const std::vector<data::CellData>& model,
const RestartConfig& restartConfig,
const int reportStepNum) {
return model;
}
}
template<class Model>
inline void
BlackoilOutputWriter::
writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellState& localWellState,
const Model& physicalModel,
bool substep)
{
// VTK output (is parallel if grid is parallel)
if( vtkWriter_ ) {
vtkWriter_->writeTimeStep( timer, localState, localWellState, false );
}
bool isIORank = output_ ;
if( parallelOutput_ && parallelOutput_->isParallel() )
{
// If this is not the initial write and no substep, then the well
// state used in the computation is actually the one of the last
// step. We need that well state for the gathering. Otherwise
// It an exception with a message like "global state does not
// contain well ..." might be thrown.
int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
(timer.reportStepNum() - 1) : timer.reportStepNum();
// collect all solutions to I/O rank
isIORank = parallelOutput_->collectToIORank( localState, localWellState, wellStateStepNumber );
}
const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
const RestartConfig& restartConfig = eclipseState_->getRestartConfig();
const int reportStepNum = timer.reportStepNum();
std::vector<data::CellData> cellData = detail::getCellData( physicalModel, restartConfig, reportStepNum );
// serial output is only done on I/O rank
if( isIORank )
{
if( asyncOutput_ ) {
// dispatch the write call to the extra thread
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, cellData, substep ) );
}
else {
// just write the data to disk
writeTimeStepSerial( timer, state, wellState, cellData, substep );
}
}
}
}
#endif

View File

@ -61,6 +61,7 @@ namespace Opm {
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
@ -120,6 +121,16 @@ namespace Opm {
private:
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (pu.phase_used[phase]) {
const int pos = pu.phase_pos[phase];
return rq_[pos].b;
}
else {
return ADB::null();
}
}
struct ReservoirResidualQuant {
ReservoirResidualQuant();