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 #1391 from totto82/RFIP

Browse Source 
Use FIP, EGRID and INIT output code in ebos.
This commit is contained in:
Atgeirr Flø Rasmussen 2018-02-06 12:54:01 +01:00 committed by GitHub
commit 532403c5fb
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 41 additions and 795 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
compareECLFiles.cmake
View File

@ -119,7 +119,7 @@ function(add_test_compare_parallel_restarted_simulation)
set(multiValueArgs TEST_ARGS) set(multiValueArgs TEST_ARGS)
cmake_parse_arguments(PARAM "$" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} ) cmake_parse_arguments(PARAM "$" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} )
set(RESULT_PATH ${BASE_RESULT_PATH}/restart/${PARAM_SIMULATOR}+${PARAM_CASENAME}) set(RESULT_PATH ${BASE_RESULT_PATH}/parallelRestart/${PARAM_SIMULATOR}+${PARAM_CASENAME})
set(TEST_ARGS ${OPM_DATA_ROOT}/${PARAM_CASENAME}/${PARAM_FILENAME} ${PARAM_TEST_ARGS}) set(TEST_ARGS ${OPM_DATA_ROOT}/${PARAM_CASENAME}/${PARAM_FILENAME} ${PARAM_TEST_ARGS})
opm_add_test(compareParallelRestartedSim_${PARAM_SIMULATOR}+${PARAM_FILENAME} NO_COMPILE opm_add_test(compareParallelRestartedSim_${PARAM_SIMULATOR}+${PARAM_FILENAME} NO_COMPILE

157
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -117,8 +117,6 @@ namespace Opm {
typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType; typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType;
//typedef typename SolutionVector :: value_type PrimaryVariables ; //typedef typename SolutionVector :: value_type PrimaryVariables ;
typedef Opm::FIPData FIPDataType;
// --------- Public methods --------- // --------- Public methods ---------
/// Construct the model. It will retain references to the /// Construct the model. It will retain references to the
@ -999,144 +997,16 @@ namespace Opm {
return computeFluidInPlace(fipnum); return computeFluidInPlace(fipnum);
} }
/// Should not be called
std::vector<std::vector<double> > std::vector<std::vector<double> >
computeFluidInPlace(const std::vector<int>& fipnum) const computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
{ {
const auto& comm = grid_.comm(); //assert(true)
const auto& gridView = ebosSimulator().gridView(); //return an empty vector
const int nc = gridView.size(/*codim=*/0); std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
int ntFip = *std::max_element(fipnum.begin(), fipnum.end());
ntFip = comm.max(ntFip);
std::vector<double> tpv(ntFip, 0.0);
std::vector<double> hcpv(ntFip, 0.0);
std::vector<std::vector<double> > regionValues(ntFip, std::vector<double>(FIPDataType::fipValues,0.0));
for (int i = 0; i<FIPDataType::fipValues; i++) {
fip_.fip[i].resize(nc,0.0);
}
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
elemIt != elemEndIt;
++elemIt)
{
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const int regionIdx = fipnum[cellIdx] - 1;
if (regionIdx < 0) {
// the given cell is not attributed to any region
continue;
}
// calculate the pore volume of the current cell. Note that the porosity
// returned by the intensive quantities is defined as the ratio of pore
// space to total cell volume and includes all pressure dependent (->
// rock compressibility) and static modifiers (MULTPV, MULTREGP, NTG,
// PORV, MINPV and friends). Also note that because of this, the porosity
// returned by the intensive quantities can be outside of the physical
// range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
}
const double b = fs.invB(phase).value();
const double s = fs.saturation(phase).value();
const unsigned flowCanonicalPhaseIdx = ebosPhaseToFlowCanonicalPhaseIdx(phase);
fip_.fip[flowCanonicalPhaseIdx][cellIdx] = b * s * pv;
regionValues[regionIdx][flowCanonicalPhaseIdx] += fip_.fip[flowCanonicalPhaseIdx][cellIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Account for gas dissolved in oil and vaporized oil
fip_.fip[FIPDataType::FIP_DISSOLVED_GAS][cellIdx] = fs.Rs().value() * fip_.fip[FIPDataType::FIP_LIQUID][cellIdx];
fip_.fip[FIPDataType::FIP_VAPORIZED_OIL][cellIdx] = fs.Rv().value() * fip_.fip[FIPDataType::FIP_VAPOUR][cellIdx];
regionValues[regionIdx][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][cellIdx];
regionValues[regionIdx][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][cellIdx];
}
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
tpv[regionIdx] += pv;
hcpv[regionIdx] += pv * hydrocarbon;
}
// sum tpv (-> total pore volume of the regions) and hcpv (-> pore volume of the
// the regions that is occupied by hydrocarbons)
comm.sum(tpv.data(), tpv.size());
comm.sum(hcpv.data(), hcpv.size());
for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const int regionIdx = fipnum[cellIdx] - 1;
if (regionIdx < 0) {
// the cell is not attributed to any region. ignore it!
continue;
}
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// calculate the pore volume of the current cell. Note that the
// porosity returned by the intensive quantities is defined as the
// ratio of pore space to total cell volume and includes all pressure
// dependent (-> rock compressibility) and static modifiers (MULTPV,
// MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
// this, the porosity returned by the intensive quantities can be
// outside of the physical range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
fip_.fip[FIPDataType::FIP_PV][cellIdx] = pv;
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
//Compute hydrocarbon pore volume weighted average pressure.
//If we have no hydrocarbon in region, use pore volume weighted average pressure instead
if (hcpv[regionIdx] > 1e-10) {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[regionIdx];
} else {
fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx] = pv * fs.pressure(FluidSystem::oilPhaseIdx).value() / tpv[regionIdx];
}
regionValues[regionIdx][FIPDataType::FIP_PV] += fip_.fip[FIPDataType::FIP_PV][cellIdx];
regionValues[regionIdx][FIPDataType::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPDataType::FIP_WEIGHTED_PRESSURE][cellIdx];
}
// sum the results over all processes
for(int regionIdx=0; regionIdx < ntFip; ++regionIdx) {
comm.sum(regionValues[regionIdx].data(), regionValues[regionIdx].size());
}
return regionValues; return regionValues;
} }
const FIPDataType& getFIPData() const {
return fip_;
}
const Simulator& ebosSimulator() const const Simulator& ebosSimulator() const
{ return ebosSimulator_; } { return ebosSimulator_; }
@ -1176,7 +1046,6 @@ namespace Opm {
std::vector<std::vector<double>> residual_norms_history_; std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_; double current_relaxation_;
BVector dx_old_; BVector dx_old_;
mutable FIPDataType fip_;
public: public:
/// return the StandardWells object /// return the StandardWells object
@ -1186,20 +1055,6 @@ namespace Opm {
const BlackoilWellModel<TypeTag>& const BlackoilWellModel<TypeTag>&
wellModel() const { return well_model_; } wellModel() const { return well_model_; }
int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const
{
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx)
return Water;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx)
return Oil;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx)
return Gas;
assert(phaseIdx < 3);
// for other phases return the index
return phaseIdx;
}
void beginReportStep() void beginReportStep()
{ {
ebosSimulator_.problem().beginEpisode(); ebosSimulator_.problem().beginEpisode();

147
opm/autodiff/BlackoilOutputEbos.hpp
View File

@ -136,13 +136,8 @@ namespace Opm
const double nextstep = -1.0, const double nextstep = -1.0,
const SimulatorReport& simulatorReport = SimulatorReport()) const SimulatorReport& simulatorReport = SimulatorReport())
{ {
data::Solution fip{};
if( output_ ) if( output_ )
{ {
// Get FIP dat
getSummaryData( fip, phaseUsage_, physicalModel, summaryConfig() );
// Add TCPU if simulatorReport is not defaulted. // Add TCPU if simulatorReport is not defaulted.
const double totalSolverTime = simulatorReport.solver_time; const double totalSolverTime = simulatorReport.solver_time;
@ -169,7 +164,7 @@ namespace Opm
const WellStateFullyImplicitBlackoil& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; const WellStateFullyImplicitBlackoil& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
// The writeOutput expects a local data::solution vector and a global data::well vector. // The writeOutput expects a local data::solution vector and a global data::well vector.
ebosSimulator_.problem().writeOutput( wellState.report(phaseUsage_), timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep, fip); ebosSimulator_.problem().writeOutput( wellState.report(phaseUsage_), timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep);
} }
} }
@ -222,157 +217,19 @@ namespace Opm
} }
} }
bool requireFIPNUM() const
{ return summaryConfig().requireFIPNUM(); }
const Grid& grid() const Grid& grid()
{ return ebosSimulator_.gridManager().grid(); } { return ebosSimulator_.gridManager().grid(); }
const Schedule& schedule() const const Schedule& schedule() const
{ return ebosSimulator_.gridManager().schedule(); } { return ebosSimulator_.gridManager().schedule(); }
const SummaryConfig& summaryConfig() const
{ return ebosSimulator_.gridManager().summaryConfig(); }
const EclipseState& eclState() const const EclipseState& eclState() const
{ return ebosSimulator_.gridManager().eclState(); } { return ebosSimulator_.gridManager().eclState(); }
bool isRestart() const { bool isRestart() const {
const auto& initconfig = eclState().getInitConfig(); const auto& initconfig = eclState().getInitConfig();
return initconfig.restartRequested(); return initconfig.restartRequested();
} }
private:
/**
* Checks if the summaryConfig has a keyword with the standardized field, region, or block prefixes.
*/
inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) {
std::string field_kw = "F" + keyword;
std::string region_kw = "R" + keyword;
std::string block_kw = "B" + keyword;
return summaryConfig.hasKeyword(field_kw)
|| summaryConfig.hasKeyword(region_kw)
|| summaryConfig.hasKeyword(block_kw);
}
/**
* Returns the data as asked for in the summaryConfig
*/
template<class Model>
void getSummaryData(data::Solution& output,
const Opm::PhaseUsage& phaseUsage,
const Model& physicalModel,
const SummaryConfig& summaryConfig) {
typedef typename Model::FIPDataType FIPDataType;
typedef typename FIPDataType::VectorType VectorType;
FIPDataType fd = physicalModel.getFIPData();
//Get shorthands for water, oil, gas
const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
/**
* Now process all of the summary config files
*/
// Water in place
if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) {
output.insert("WIP",
Opm::UnitSystem::measure::volume,
std::move( fd.fip[ FIPDataType::FIP_AQUA ] ),
data::TargetType::SUMMARY );
}
if (liquid_active) {
const VectorType& oipl = fd.fip[FIPDataType::FIP_LIQUID];
VectorType oip ( oipl );
const size_t size = oip.size();
const VectorType& oipg = vapour_active ? fd.fip[FIPDataType::FIP_VAPORIZED_OIL] : VectorType(size, 0.0);
if( vapour_active )
{
// oip = oipl + oipg
for( size_t i=0; i<size; ++ i ) {
oip[ i ] += oipg[ i ];
}
}
//Oil in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "OIPL")) {
output.insert("OIPL",
Opm::UnitSystem::measure::volume,
std::move( oipl ),
data::TargetType::SUMMARY );
}
//Oil in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "OIPG")) {
output.insert("OIPG",
Opm::UnitSystem::measure::volume,
std::move( oipg ),
data::TargetType::SUMMARY );
}
// Oil in place (in liquid and gas phases)
if (hasFRBKeyword(summaryConfig, "OIP") || hasFRBKeyword(summaryConfig, "OE")) {
output.insert("OIP",
Opm::UnitSystem::measure::volume,
std::move( oip ),
data::TargetType::SUMMARY );
}
}
if (vapour_active) {
const VectorType& gipg = fd.fip[ FIPDataType::FIP_VAPOUR];
VectorType gip( gipg );
const size_t size = gip.size();
const VectorType& gipl = liquid_active ? fd.fip[ FIPDataType::FIP_DISSOLVED_GAS ] : VectorType(size,0.0);
if( liquid_active )
{
// gip = gipg + gipl
for( size_t i=0; i<size; ++ i ) {
gip[ i ] += gipl[ i ];
}
}
// Gas in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "GIPG")) {
output.insert("GIPG",
Opm::UnitSystem::measure::volume,
std::move( gipg ),
data::TargetType::SUMMARY );
}
// Gas in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "GIPL")) {
output.insert("GIPL",
Opm::UnitSystem::measure::volume,
std::move( gipl ),
data::TargetType::SUMMARY );
}
// Gas in place (in both liquid and gas phases)
if (hasFRBKeyword(summaryConfig, "GIP")) {
output.insert("GIP",
Opm::UnitSystem::measure::volume,
std::move( gip ),
data::TargetType::SUMMARY );
}
}
// Cell pore volume in reservoir conditions
if (hasFRBKeyword(summaryConfig, "RPV")) {
output.insert("RPV",
Opm::UnitSystem::measure::volume,
std::move( fd.fip[FIPDataType::FIP_PV]),
data::TargetType::SUMMARY );
}
// Pressure averaged value (hydrocarbon pore volume weighted)
if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) {
output.insert("PRH",
Opm::UnitSystem::measure::pressure,
std::move(fd.fip[FIPDataType::FIP_WEIGHTED_PRESSURE]),
data::TargetType::SUMMARY );
}
}
protected: protected:
const bool output_; const bool output_;

294
opm/autodiff/FlowMainEbos.hpp
View File

@ -26,9 +26,7 @@
#include <sys/utsname.h> #include <sys/utsname.h>
#include <opm/simulators/ParallelFileMerger.hpp> #include <opm/simulators/ParallelFileMerger.hpp>
#include <opm/simulators/ensureDirectoryExists.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp> #include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp> #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
@ -110,7 +108,6 @@ namespace Opm
printPRTHeader(); printPRTHeader();
extractMessages(); extractMessages();
runDiagnostics(); runDiagnostics();
writeInit();
setupOutputWriter(); setupOutputWriter();
setupLinearSolver(); setupLinearSolver();
createSimulator(); createSimulator();
@ -257,23 +254,7 @@ namespace Opm
<< std::endl; << std::endl;
} }
output_to_files_ = output_cout_ && output_ > OUTPUT_NONE; output_to_files_ = output_cout_ && (output_ != OUTPUT_NONE);
// Setup output directory.
auto& ioConfig = eclState().getIOConfig();
// Default output directory is the directory where the deck is found.
const std::string default_output_dir = ioConfig.getOutputDir();
output_dir_ = param_.getDefault("output_dir", default_output_dir);
// Override output directory if user specified.
ioConfig.setOutputDir(output_dir_);
// Write parameters used for later reference. (only if rank is zero)
if (output_to_files_) {
// Create output directory if needed.
ensureDirectoryExists(output_dir_);
// Write simulation parameters.
param_.writeParam(output_dir_ + "/simulation.param");
}
} }
// Setup OpmLog backend with output_dir. // Setup OpmLog backend with output_dir.
@ -413,9 +394,22 @@ namespace Opm
argv.push_back("flow_ebos"); argv.push_back("flow_ebos");
std::string deckFileParam("--ecl-deck-file-name="); std::string deckFileParam("--ecl-deck-file-name=");
deckFileParam += param_.get<std::string>("deck_filename"); const std::string& deckFileName = param_.get<std::string>("deck_filename");
deckFileParam += deckFileName;
argv.push_back(deckFileParam.c_str()); argv.push_back(deckFileParam.c_str());
const std::string default_output_dir = boost::filesystem::basename(deckFileName);
output_dir_ = param_.getDefault("output_dir", default_output_dir);
std::string outputDirParam("--ecl-output-dir=");
outputDirParam += output_dir_;
argv.push_back(outputDirParam.c_str());
const bool restart_double_si = param_.getDefault("restart_double_si", false);
std::string outputDoublePrecisionParam("--ecl-output-double-precision=");
outputDoublePrecisionParam += restart_double_si ? "true" : "false";
argv.push_back(outputDoublePrecisionParam.c_str());
#if defined(_OPENMP) #if defined(_OPENMP)
std::string numThreadsParam("--threads-per-process="); std::string numThreadsParam("--threads-per-process=");
int numThreads = omp_get_max_threads(); int numThreads = omp_get_max_threads();
@ -430,10 +424,6 @@ namespace Opm
ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false)); ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false));
ebosSimulator_->model().applyInitialSolution(); ebosSimulator_->model().applyInitialSolution();
// Create a grid with a global view.
globalGrid_.reset(new Grid(grid()));
globalGrid_->switchToGlobalView();
try { try {
if (output_cout_) { if (output_cout_) {
MissingFeatures::checkKeywords(deck()); MissingFeatures::checkKeywords(deck());
@ -474,9 +464,6 @@ namespace Opm
const Schedule& schedule() const const Schedule& schedule() const
{ return ebosSimulator_->gridManager().schedule(); } { return ebosSimulator_->gridManager().schedule(); }
const SummaryConfig& summaryConfig() const
{ return ebosSimulator_->gridManager().summaryConfig(); }
// Extract messages from parser. // Extract messages from parser.
// Writes to: // Writes to:
// OpmLog singleton. // OpmLog singleton.
@ -523,27 +510,6 @@ namespace Opm
diagnostic.diagnosis(eclState(), deck(), this->grid()); diagnostic.diagnosis(eclState(), deck(), this->grid());
} }
void writeInit()
{
bool output = ( output_ > OUTPUT_LOG_ONLY );
bool output_ecl = param_.getDefault("output_ecl", true);
auto int_vectors = computeCellRanks(output, output_ecl);
if( output && output_ecl && grid().comm().rank() == 0 )
{
exportNncStructure_();
const EclipseGrid& inputGrid = eclState().getInputGrid();
eclIO_.reset(new EclipseIO(eclState(),
UgGridHelpers::createEclipseGrid( this->globalGrid() , inputGrid ),
schedule(),
summaryConfig()));
eclIO_->writeInitial(computeLegacySimProps_(), int_vectors, nnc_);
Problem& problem = ebosProblem();
problem.setEclIO(std::move(eclIO_));
}
}
// Setup output writer. // Setup output writer.
// Writes to: // Writes to:
// output_writer_ // output_writer_
@ -689,232 +655,6 @@ namespace Opm
Grid& grid() Grid& grid()
{ return ebosSimulator_->gridManager().grid(); } { return ebosSimulator_->gridManager().grid(); }
const Grid& globalGrid()
{ return *globalGrid_; }
Problem& ebosProblem()
{ return ebosSimulator_->problem(); }
const Problem& ebosProblem() const
{ return ebosSimulator_->problem(); }
std::shared_ptr<MaterialLawManager> materialLawManager()
{ return ebosProblem().materialLawManager(); }
Scalar gravity() const
{ return ebosProblem().gravity()[2]; }
std::map<std::string, std::vector<int> > computeCellRanks(bool output, bool output_ecl)
{
std::map<std::string, std::vector<int> > integerVectors;
if( output && output_ecl && grid().comm().size() > 1 )
{
typedef typename Grid::LeafGridView GridView;
#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
#else
// Get the owner rank number for each cell
using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout>;
#endif
using Handle = CellOwnerDataHandle<ElementMapper>;
const Grid& globalGrid = this->globalGrid();
const auto& globalGridView = globalGrid.leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
ElementMapper globalMapper(globalGridView, Dune::mcmgElementLayout());
#else
ElementMapper globalMapper(globalGridView);
#endif
const auto globalSize = globalGrid.size(0);
std::vector<int> ranks(globalSize, -1);
Handle handle(globalMapper, ranks);
this->grid().gatherData(handle);
integerVectors.emplace("MPI_RANK", ranks);
}
return integerVectors;
}
data::Solution computeLegacySimProps_()
{
const int* dims = UgGridHelpers::cartDims(grid());
const int globalSize = dims[0]*dims[1]*dims[2];
data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
for (size_t i = 0; i < tranx.data.size(); ++i) {
tranx.data[0] = 0.0;
trany.data[0] = 0.0;
tranz.data[0] = 0.0;
}
const Grid& globalGrid = this->globalGrid();
const auto& globalGridView = globalGrid.leafGridView();
typedef typename Grid::LeafGridView GridView;
#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
#endif
const auto& cartesianCellIdx = globalGrid.globalCell();
const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
if (grid().comm().size() < 2) {
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
}
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
if (!is.neighbor())
{
continue; // intersection is on the domain boundary
}
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
{
continue; // we only need to handle each connection once, thank you.
}
int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]);
int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]);
if (gc2 - gc1 == 1) {
tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
}
if (gc2 - gc1 == dims[0]) {
trany.data[gc1] = globalTrans->transmissibility(c1, c2);
}
if (gc2 - gc1 == dims[0]*dims[1]) {
tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
}
}
}
return {{"TRANX" , tranx},
{"TRANY" , trany} ,
{"TRANZ" , tranz}};
}
void exportNncStructure_()
{
nnc_ = eclState().getInputNNC();
int nx = eclState().getInputGrid().getNX();
int ny = eclState().getInputGrid().getNY();
//int nz = eclState().getInputGrid().getNZ()
const Grid& globalGrid = this->globalGrid();
const auto& globalGridView = globalGrid.leafGridView();
typedef typename Grid::LeafGridView GridView;
#if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
#endif
const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
if (grid().comm().size() < 2) {
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
}
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
if (!is.neighbor())
{
continue; // intersection is on the domain boundary
}
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
{
continue; // we only need to handle each connection once, thank you.
}
// TODO (?): use the cartesian index mapper to make this code work
// with grids other than Dune::CpGrid. The problem is that we need
// the a mapper for the sequential grid, not for the distributed one.
int cc1 = globalGrid.globalCell()[c1];
int cc2 = globalGrid.globalCell()[c2];
if (std::abs(cc1 - cc2) != 1 &&
std::abs(cc1 - cc2) != nx &&
std::abs(cc1 - cc2) != nx*ny)
{
nnc_.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2));
}
}
}
}
/// Convert saturations from a vector of individual phase saturation vectors
/// to an interleaved format where all values for a given cell come before all
/// values for the next cell, all in a single vector.
template <class FluidSystem>
void convertSats(std::vector<double>& sat_interleaved, const std::vector< std::vector<double> >& sat, const PhaseUsage& pu)
{
assert(sat.size() == 3);
const auto nc = sat[0].size();
const auto np = sat_interleaved.size() / nc;
for (size_t c = 0; c < nc; ++c) {
if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sat_p = sat[ FluidSystem::oilPhaseIdx];
sat_interleaved[np*c + opos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const std::vector<double>& sat_p = sat[ FluidSystem::waterPhaseIdx];
sat_interleaved[np*c + wpos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const std::vector<double>& sat_p = sat[ FluidSystem::gasPhaseIdx];
sat_interleaved[np*c + gpos] = sat_p[c];
}
}
}
std::unique_ptr<EbosSimulator> ebosSimulator_; std::unique_ptr<EbosSimulator> ebosSimulator_;
int mpi_rank_ = 0; int mpi_rank_ = 0;
bool output_cout_ = false; bool output_cout_ = false;
@ -923,15 +663,11 @@ namespace Opm
ParameterGroup param_; ParameterGroup param_;
bool output_to_files_ = false; bool output_to_files_ = false;
std::string output_dir_ = std::string("."); std::string output_dir_ = std::string(".");
NNC nnc_;
std::unique_ptr<EclipseIO> eclIO_;
std::unique_ptr<OutputWriter> output_writer_; std::unique_ptr<OutputWriter> output_writer_;
boost::any parallel_information_; boost::any parallel_information_;
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_; std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
std::unique_ptr<Simulator> simulator_; std::unique_ptr<Simulator> simulator_;
std::string logFile_; std::string logFile_;
// Needs to be shared pointer because it gets initialzed before MPI_Init.
std::shared_ptr<Grid> globalGrid_;
}; };
} // namespace Opm } // namespace Opm

227
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -31,7 +31,6 @@
#include <opm/autodiff/BlackoilWellModel.hpp> #include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/moduleVersion.hpp> #include <opm/autodiff/moduleVersion.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp> #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/core/utility/initHydroCarbonState.hpp>
#include <opm/core/utility/StopWatch.hpp> #include <opm/core/utility/StopWatch.hpp>
#include <opm/common/Exceptions.hpp> #include <opm/common/Exceptions.hpp>
@ -119,7 +118,6 @@ public:
is_parallel_run_ = ( info.communicator().size() > 1 ); is_parallel_run_ = ( info.communicator().size() > 1 );
} }
#endif #endif
createLocalFipnum();
} }
/// Run the simulation. /// Run the simulation.
@ -198,20 +196,11 @@ public:
auto solver = createSolver(well_model); auto solver = createSolver(well_model);
// Compute orignal fluid in place if this has not been done yet
if (originalFluidInPlace_.data.empty()) {
originalFluidInPlace_ = computeFluidInPlace(*solver);
}
// write the inital state at the report stage // write the inital state at the report stage
if (timer.initialStep()) { if (timer.initialStep()) {
Dune::Timer perfTimer; Dune::Timer perfTimer;
perfTimer.start(); perfTimer.start();
if (terminal_output_) {
outputFluidInPlace(timer, originalFluidInPlace_);
}
// No per cell data is written for initial step, but will be // No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating // for subsequent steps, when we have started simulating
output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() ); output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() );
@ -250,8 +239,7 @@ public:
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum()); events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, nullptr );
output_writer_.requireFIPNUM() ? &fipnum_ : nullptr );
report += stepReport; report += stepReport;
failureReport_ += adaptiveTimeStepping->failureReport(); failureReport_ += adaptiveTimeStepping->failureReport();
} }
@ -288,26 +276,30 @@ public:
// Increment timer, remember well state. // Increment timer, remember well state.
++timer; ++timer;
// Compute current fluid in place.
const auto currentFluidInPlace = computeFluidInPlace(*solver);
if (terminal_output_ ) if (terminal_output_ )
{ {
outputFluidInPlace(timer, currentFluidInPlace); if (!timer.initialStep()) {
const std::string version = moduleVersionName();
std::string msg = outputTimestampFIP(timer, version);
"Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; " }
"total solver time " + std::to_string(report.solver_time) + " seconds.";
OpmLog::debug(msg);
} }
// write simulation state at the report stage // write simulation state at the report stage
Dune::Timer perfTimer; Dune::Timer perfTimer;
perfTimer.start(); perfTimer.start();
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0; const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report); output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report);
report.output_write_time += perfTimer.stop(); report.output_write_time += perfTimer.stop();
if (terminal_output_ )
{
std::string msg =
"Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; "
"total solver time " + std::to_string(report.solver_time) + " seconds.";
OpmLog::debug(msg);
}
} }
// Stop timer and create timing report // Stop timer and create timing report
@ -337,150 +329,6 @@ protected:
return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model))); return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
} }
void createLocalFipnum()
{
const std::vector<int>& fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
// Get compressed cell fipnum.
fipnum_.resize(Opm::UgGridHelpers::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[Opm::UgGridHelpers::globalCell(grid())[c]];
}
}
}
void FIPUnitConvert(const UnitSystem& units,
std::vector<std::vector<double>>& fip)
{
for (size_t i = 0; i < fip.size(); ++i) {
FIPUnitConvert(units, fip[i]);
}
}
void FIPUnitConvert(const UnitSystem& units,
std::vector<double>& fip)
{
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
fip[0] = unit::convert::to(fip[0], unit::stb);
fip[1] = unit::convert::to(fip[1], unit::stb);
fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
fip[4] = unit::convert::to(fip[4], unit::stb);
fip[5] = unit::convert::to(fip[5], unit::stb);
fip[6] = unit::convert::to(fip[6], unit::psia);
}
else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
fip[6] = unit::convert::to(fip[6], unit::barsa);
}
else {
OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
}
}
std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip)
{
std::vector<double> totals(7,0.0);
for (int i = 0; i < 5; ++i) {
for (size_t reg = 0; reg < fip.size(); ++reg) {
totals[i] += fip[reg][i];
}
}
const auto& gridView = ebosSimulator_.gridManager().gridView();
const auto& comm = gridView.comm();
double pv_hydrocarbon_sum = 0.0;
double p_pv_hydrocarbon_sum = 0.0;
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
// calculate the pore volume of the current cell. Note that the
// porosity returned by the intensive quantities is defined as the
// ratio of pore space to total cell volume and includes all pressure
// dependent (-> rock compressibility) and static modifiers (MULTPV,
// MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
// this, the porosity returned by the intensive quantities can be
// outside of the physical range [0, 1] in pathetic cases.
const double pv =
ebosSimulator_.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
totals[5] += pv;
pv_hydrocarbon_sum += pv*hydrocarbon;
p_pv_hydrocarbon_sum += p*pv*hydrocarbon;
}
pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum);
p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum);
totals[5] = comm.sum(totals[5]);
totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum);
return totals;
}
struct FluidInPlace
{
std::vector<std::vector<double>> data;
std::vector<double> totals;
};
FluidInPlace computeFluidInPlace(const Solver& solver)
{
FluidInPlace fip;
fip.data = solver.computeFluidInPlace(fipnum_);
fip.totals = FIPTotals(fip.data);
FIPUnitConvert(eclState().getUnits(), fip.data);
FIPUnitConvert(eclState().getUnits(), fip.totals);
return fip;
}
void outputFluidInPlace(const SimulatorTimer& timer,
const FluidInPlace& currentFluidInPlace)
{
if (!timer.initialStep()) {
const std::string version = moduleVersionName();
outputTimestampFIP(timer, version);
}
outputRegionFluidInPlace(originalFluidInPlace_.totals,
currentFluidInPlace.totals,
eclState().getUnits(),
0);
for (size_t reg = 0; reg < originalFluidInPlace_.data.size(); ++reg) {
outputRegionFluidInPlace(originalFluidInPlace_.data[reg],
currentFluidInPlace.data[reg],
eclState().getUnits(),
reg+1);
}
}
void outputTimestampFIP(const SimulatorTimer& timer, const std::string version) void outputTimestampFIP(const SimulatorTimer& timer, const std::string version)
{ {
std::ostringstream ss; std::ostringstream ss;
@ -495,50 +343,6 @@ protected:
OpmLog::note(ss.str()); OpmLog::note(ss.str());
} }
void outputRegionFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
{
std::ostringstream ss;
if (!reg) {
ss << " ===================================================\n"
<< " : Field Totals :\n";
} else {
ss << " ===================================================\n"
<< " : FIPNUM report region "
<< std::setw(2) << reg << " :\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Porv volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
}
ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
<< ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
<< std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
<< ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
<< std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
<< ":========================:==========================================:================:==========================================:\n";
OpmLog::note(ss.str());
}
const EclipseState& eclState() const const EclipseState& eclState() const
{ return ebosSimulator_.gridManager().eclState(); } { return ebosSimulator_.gridManager().eclState(); }
@ -550,9 +354,6 @@ protected:
// Data. // Data.
Simulator& ebosSimulator_; Simulator& ebosSimulator_;
std::vector<int> fipnum_;
FluidInPlace originalFluidInPlace_;
typedef typename Solver::SolverParameters SolverParameters; typedef typename Solver::SolverParameters SolverParameters;
SimulatorReport failureReport_; SimulatorReport failureReport_;

7
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
View File

@ -333,11 +333,6 @@ namespace Opm
if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) { if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
} else { } else {
if ( timer.initialStep() )
{
// Set the initial OIP
eclIO_->overwriteInitialOIP(simProps);
}
// ... insert "extra" data (KR, VISC, ...) // ... insert "extra" data (KR, VISC, ...)
const int reportStepForOutput = substep ? timer.reportStepNum() + 1 : timer.reportStepNum(); const int reportStepForOutput = substep ? timer.reportStepNum() + 1 : timer.reportStepNum();
eclIO_->writeTimeStep(reportStepForOutput, eclIO_->writeTimeStep(reportStepForOutput,
@ -346,6 +341,8 @@ namespace Opm
simProps, simProps,
wellState.report(phaseUsage_), wellState.report(phaseUsage_),
miscSummaryData, miscSummaryData,
{}, //regionData
{}, //blockData
extraRestartData, extraRestartData,
restart_double_si_); restart_double_si_);
} }

2
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
View File

@ -333,7 +333,7 @@ namespace Opm {
perfTimer.start(); perfTimer.start();
bool substep = true; bool substep = true;
const auto& physicalModel = solver.model(); const auto& physicalModel = solver.model();
outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep); outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep, -1.0, substepReport);
report.output_write_time += perfTimer.secsSinceStart(); report.output_write_time += perfTimer.secsSinceStart();
} }