diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0bfcff618..6d0afdb6c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -210,6 +210,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellDensitySegmented.hpp opm/autodiff/WellStateFullyImplicitBlackoil.hpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp + opm/autodiff/BlackoilOutputEbos.hpp opm/autodiff/VFPProperties.hpp opm/autodiff/VFPHelpers.hpp opm/autodiff/VFPProdProperties.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 567e2db9c..d331af89d 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -1125,315 +1125,6 @@ namespace Opm { return regionValues; } - SimulationDataContainer getSimulatorData ( const SimulationDataContainer& /*localState*/) const - { - typedef std::vector VectorType; - - const auto& ebosModel = ebosSimulator().model(); - const auto& phaseUsage = phaseUsage_; - - // extract everything which can possibly be written to disk - const int numCells = ebosModel.numGridDof(); - const int num_phases = numPhases(); - - SimulationDataContainer simData( numCells, 0, num_phases ); - - //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]; - - const int aqua_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Aqua ]; - const int liquid_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Liquid ]; - const int vapour_pos = phaseUsage.phase_pos[ Opm::PhaseUsage::Vapour ]; - - VectorType zero; - - VectorType& pressureOil = simData.pressure(); - VectorType& temperature = simData.temperature(); - VectorType& saturation = simData.saturation(); - - // WATER - if( aqua_active ) { - simData.registerCellData( "1OVERBW", 1 ); - simData.registerCellData( "WAT_DEN", 1 ); - simData.registerCellData( "WAT_VISC", 1 ); - simData.registerCellData( "WATKR", 1 ); - } - - VectorType& bWater = aqua_active ? simData.getCellData( "1OVERBW" ) : zero; - VectorType& rhoWater = aqua_active ? simData.getCellData( "WAT_DEN" ) : zero; - VectorType& muWater = aqua_active ? simData.getCellData( "WAT_VISC" ) : zero; - VectorType& krWater = aqua_active ? simData.getCellData( "WATKR" ) : zero; - - // OIL - if( liquid_active ) { - simData.registerCellData( "1OVERBO", 1 ); - simData.registerCellData( "OIL_DEN", 1 ); - simData.registerCellData( "OIL_VISC", 1 ); - simData.registerCellData( "OILKR", 1 ); - } - - VectorType& bOil = liquid_active ? simData.getCellData( "1OVERBO" ) : zero; - VectorType& rhoOil = liquid_active ? simData.getCellData( "OIL_DEN" ) : zero; - VectorType& muOil = liquid_active ? simData.getCellData( "OIL_VISC" ) : zero; - VectorType& krOil = liquid_active ? simData.getCellData( "OILKR" ) : zero; - - // GAS - if( vapour_active ) { - simData.registerCellData( "1OVERBG", 1 ); - simData.registerCellData( "GAS_DEN", 1 ); - simData.registerCellData( "GAS_VISC", 1 ); - simData.registerCellData( "GASKR", 1 ); - } - - VectorType& bGas = vapour_active ? simData.getCellData( "1OVERBG" ) : zero; - VectorType& rhoGas = vapour_active ? simData.getCellData( "GAS_DEN" ) : zero; - VectorType& muGas = vapour_active ? simData.getCellData( "GAS_VISC" ) : zero; - VectorType& krGas = vapour_active ? simData.getCellData( "GASKR" ) : zero; - - simData.registerCellData( BlackoilState::GASOILRATIO, 1 ); - simData.registerCellData( BlackoilState::RV, 1 ); - simData.registerCellData( "RSSAT", 1 ); - simData.registerCellData( "RVSAT", 1 ); - - VectorType& Rs = simData.getCellData( BlackoilState::GASOILRATIO ); - VectorType& Rv = simData.getCellData( BlackoilState::RV ); - VectorType& RsSat = simData.getCellData( "RSSAT" ); - VectorType& RvSat = simData.getCellData( "RVSAT" ); - - simData.registerCellData( "PBUB", 1 ); - simData.registerCellData( "PDEW", 1 ); - - VectorType& Pb = simData.getCellData( "PBUB" ); - VectorType& Pd = simData.getCellData( "PDEW" ); - - simData.registerCellData( "SOMAX", 1 ); - VectorType& somax = simData.getCellData( "SOMAX" ); - - // Two components for hysteresis parameters - // pcSwMdc/krnSwMdc, one for oil-water and one for gas-oil - simData.registerCellData( "PCSWMDC_GO", 1 ); - simData.registerCellData( "KRNSWMDC_GO", 1 ); - - simData.registerCellData( "PCSWMDC_OW", 1 ); - simData.registerCellData( "KRNSWMDC_OW", 1 ); - - VectorType& pcSwMdc_go = simData.getCellData( "PCSWMDC_GO" ); - VectorType& krnSwMdc_go = simData.getCellData( "KRNSWMDC_GO" ); - - VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" ); - VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" ); - - if (has_solvent_) { - simData.registerCellData( "SSOL", 1 ); - } - VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero; - - if (has_polymer_) { - simData.registerCellData( "POLYMER", 1 ); - } - VectorType& cpolymer = has_polymer_ ? simData.getCellData( "POLYMER" ) : zero; - - std::vector failed_cells_pb; - std::vector failed_cells_pd; - const auto& gridView = ebosSimulator().gridView(); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - ElementContext elemCtx(ebosSimulator()); - - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - - 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 int satIdx = cellIdx * num_phases; - - pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value(); - - temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value(); - - somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx); - - const auto& matLawManager = ebosSimulator().problem().materialLawManager(); - if (matLawManager->enableHysteresis()) { - matLawManager->oilWaterHysteresisParams( - pcSwMdc_ow[cellIdx], - krnSwMdc_ow[cellIdx], - cellIdx); - matLawManager->gasOilHysteresisParams( - pcSwMdc_go[cellIdx], - krnSwMdc_go[cellIdx], - cellIdx); - } - - if (aqua_active) { - saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value(); - bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value(); - rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value(); - muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value(); - krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value(); - } - if (vapour_active) { - saturation[ satIdx + vapour_pos ] = fs.saturation(FluidSystem::gasPhaseIdx).value(); - bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value(); - rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value(); - muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value(); - krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value(); - Rs[cellIdx] = fs.Rs().value(); - Rv[cellIdx] = fs.Rv().value(); - RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, - FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex(), - /*maxOilSaturation=*/1.0).value(); - RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, - FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex(), - /*maxOilSaturation=*/1.0).value(); - try { - Pb[cellIdx] = FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()).value(); - } - catch (const NumericalProblem& e) { - const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx]; - failed_cells_pb.push_back(globalIdx); - } - try { - Pd[cellIdx] = FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()).value(); - } - catch (const NumericalProblem& e) { - const auto globalIdx = ebosSimulator_.gridManager().grid().globalCell()[cellIdx]; - failed_cells_pd.push_back(globalIdx); - } - } - if( liquid_active ) - { - saturation[ satIdx + liquid_pos ] = fs.saturation(FluidSystem::oilPhaseIdx).value(); - bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value(); - rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value(); - muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value(); - krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value(); - } - - if (has_solvent_) - { - ssol[cellIdx] = intQuants.solventSaturation().value(); - } - - if (has_polymer_) - { - cpolymer[cellIdx] = intQuants.polymerConcentration().value(); - } - - // hack to make the intial output of rs and rv Ecl compatible. - // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step - // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite - // rs and rv with the values computed in the initially. - // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values. - if (ebosSimulator_.episodeIndex() < 0 && vapour_active && liquid_active ) { - - const auto& fs_updated = ebosSimulator().problem().initialFluidState(cellIdx); - - // use initial rs and rv values - Rv[cellIdx] = fs_updated.Rv(); - Rs[cellIdx] = fs_updated.Rs(); - - //re-compute the volume factors, viscosities and densities. - rhoOil[cellIdx] = FluidSystem::density(fs_updated, - FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()); - rhoGas[cellIdx] = FluidSystem::density(fs_updated, - FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()); - - bOil[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated, - FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()); - bGas[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated, - FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()); - muOil[cellIdx] = FluidSystem::viscosity(fs_updated, - FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()); - muGas[cellIdx] = FluidSystem::viscosity(fs_updated, - FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()); - - } - } - - const size_t max_num_cells_faillog = 20; - - int pb_size = failed_cells_pb.size(), pd_size = failed_cells_pd.size(); - std::vector displ_pb, displ_pd, recv_len_pb, recv_len_pd; - const auto& comm = grid_.comm(); - - if ( comm.rank() == 0 ) - { - displ_pb.resize(comm.size()+1, 0); - displ_pd.resize(comm.size()+1, 0); - recv_len_pb.resize(comm.size()); - recv_len_pd.resize(comm.size()); - } - - comm.gather(&pb_size, recv_len_pb.data(), 1, 0); - comm.gather(&pd_size, recv_len_pd.data(), 1, 0); - std::partial_sum(recv_len_pb.begin(), recv_len_pb.end(), displ_pb.begin()+1); - std::partial_sum(recv_len_pd.begin(), recv_len_pd.end(), displ_pd.begin()+1); - std::vector global_failed_cells_pb, global_failed_cells_pd; - - if ( comm.rank() == 0 ) - { - global_failed_cells_pb.resize(displ_pb.back()); - global_failed_cells_pd.resize(displ_pd.back()); - } - - comm.gatherv(failed_cells_pb.data(), static_cast(failed_cells_pb.size()), - global_failed_cells_pb.data(), recv_len_pb.data(), - displ_pb.data(), 0); - comm.gatherv(failed_cells_pd.data(), static_cast(failed_cells_pd.size()), - global_failed_cells_pd.data(), recv_len_pd.data(), - displ_pd.data(), 0); - std::sort(global_failed_cells_pb.begin(), global_failed_cells_pb.end()); - std::sort(global_failed_cells_pd.begin(), global_failed_cells_pd.end()); - - if (global_failed_cells_pb.size() > 0) { - std::stringstream errlog; - errlog << "Finding the bubble point pressure failed for " << global_failed_cells_pb.size() << " cells ["; - errlog << global_failed_cells_pb[0]; - const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size()); - for (size_t i = 1; i < max_elems; ++i) { - errlog << ", " << global_failed_cells_pb[i]; - } - if (global_failed_cells_pb.size() > max_num_cells_faillog) { - errlog << ", ..."; - } - errlog << "]"; - OpmLog::warning("Bubble point numerical problem", errlog.str()); - } - if (global_failed_cells_pd.size() > 0) { - std::stringstream errlog; - errlog << "Finding the dew point pressure failed for " << global_failed_cells_pd.size() << " cells ["; - errlog << global_failed_cells_pd[0]; - const size_t max_elems = std::min(max_num_cells_faillog, global_failed_cells_pd.size()); - for (size_t i = 1; i < max_elems; ++i) { - errlog << ", " << global_failed_cells_pd[i]; - } - if (global_failed_cells_pd.size() > max_num_cells_faillog) { - errlog << ", ..."; - } - errlog << "]"; - OpmLog::warning("Dew point numerical problem", errlog.str()); - } - - return simData; - } - const FIPDataType& getFIPData() const { return fip_; } diff --git a/opm/autodiff/BlackoilOutputEbos.hpp b/opm/autodiff/BlackoilOutputEbos.hpp new file mode 100644 index 000000000..fccf45500 --- /dev/null +++ b/opm/autodiff/BlackoilOutputEbos.hpp @@ -0,0 +1,391 @@ +/* + Copyright (c) 2017 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_BLACKOILOUTPUTEBOS_HEADER_INCLUDED +#define OPM_BLACKOILOUTPUTEBOS_HEADER_INCLUDED + + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#ifdef HAVE_OPM_GRID +#include +#endif +namespace Opm +{ + + + /// Extra data to read/write for OPM restarting + struct ExtraData + { + double suggested_step = -1.0; + }; + + + /** \brief Wrapper ECL output. */ + template + class BlackoilOutputEbos + { + public: + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + // constructor creating different sub writers + BlackoilOutputEbos(Simulator& ebosSimulator, + const ParameterGroup& param) + : output_( [ ¶m ] () -> bool { + // If output parameter is true or all, then we do output + const std::string outputString = param.getDefault("output", std::string("all")); + return ( outputString == "all" || outputString == "true" ); + }() + ), + ebosSimulator_(ebosSimulator), + phaseUsage_(phaseUsageFromDeck(eclState())), + parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid(), eclState(), schedule(), phaseUsage_.num_phases, phaseUsage_ ) : 0 ), + restart_double_si_( output_ ? param.getDefault("restart_double_si", false) : false ), + asyncOutput_() + { + // For output. + if ( output_ ) + { + + // create output thread if enabled and rank is I/O rank + // async output is enabled by default if pthread are enabled + #if HAVE_PTHREAD + const bool asyncOutputDefault = true; + #else + const bool asyncOutputDefault = false; + #endif + if( param.getDefault("async_output", asyncOutputDefault ) ) + { + const bool isIORank = parallelOutput_ ? parallelOutput_->isIORank() : true; + #if HAVE_PTHREAD + asyncOutput_.reset( new ThreadHandle( isIORank ) ); + #else + OPM_THROW(std::runtime_error,"Pthreads were not found, cannot enable async_output"); + #endif + } + } + } + + + + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight. This function will extract the + * requested output cell properties specified by the RPTRST keyword + * and write these to file. + */ + template + void writeTimeStep(const SimulatorTimerInterface& timer, + const SimulationDataContainer& reservoirStateDummy, + const Opm::WellStateFullyImplicitBlackoil& /*wellStateDummy*/, + const Model& physicalModel, + const bool substep = false, + const double nextstep = -1.0, + const SimulatorReport& simulatorReport = SimulatorReport()) + { + data::Solution fip{}; + + if( output_ ) + { + // Get FIP dat + getSummaryData( fip, phaseUsage_, physicalModel, summaryConfig() ); + + // Add TCPU if simulatorReport is not defaulted. + const double totalSolverTime = simulatorReport.solver_time; + + const Opm::WellStateFullyImplicitBlackoil& localWellState = physicalModel.wellModel().wellState(); + + 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. + // The distribution of data::solution is not done here + data::Solution localCellDataDummy{}; + int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ? + (timer.reportStepNum() - 1) : timer.reportStepNum(); + // collect all solutions to I/O rank + parallelOutput_->collectToIORank( reservoirStateDummy, localWellState, + localCellDataDummy, + wellStateStepNumber ); + // Note that at this point the extraData are assumed to be global, i.e. identical across all processes. + } + + const WellStateFullyImplicitBlackoil& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; + + // 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); + } + } + + template + void initFromRestartFile(const PhaseUsage& /*phaseUsage*/, + const Grid& /*grid */, + SimulationDataContainer& simulatorstate, + WellState& wellstate, + ExtraData& extra) { + + std::map extra_keys { + {"OPMEXTRA" , false} + }; + + // gives a dummy dynamic_list_econ_limited + DynamicListEconLimited dummy_list_econ_limited; + const auto& defunct_well_names = ebosSimulator_.gridManager().defunctWellNames(); + WellsManager wellsmanager(eclState(), + schedule(), + eclState().getInitConfig().getRestartStep(), + Opm::UgGridHelpers::numCells(grid()), + Opm::UgGridHelpers::globalCell(grid()), + Opm::UgGridHelpers::cartDims(grid()), + Opm::UgGridHelpers::dimensions(grid()), + Opm::UgGridHelpers::cell2Faces(grid()), + Opm::UgGridHelpers::beginFaceCentroids(grid()), + dummy_list_econ_limited, + grid().comm().size() > 1, + defunct_well_names); + + const Wells* wells = wellsmanager.c_wells(); + + std::map solution_keys {}; + auto restart_values = ebosSimulator_.problem().eclIO().loadRestart(solution_keys, extra_keys); + + const int nw = wells->number_of_wells; + if (nw > 0) { + wellstate.resize(wells, simulatorstate, phaseUsage_ ); //Resize for restart step + wellsToState( restart_values.wells, phaseUsage_, wellstate ); + } + + const auto opmextra_iter = restart_values.extra.find("OPMEXTRA"); + if (opmextra_iter != restart_values.extra.end()) { + std::vector opmextra = opmextra_iter->second; + assert(opmextra.size() == 1); + extra.suggested_step = opmextra[0]; + } else { + OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run."); + extra.suggested_step = -1.0; + } + } + + bool requireFIPNUM() const + { return summaryConfig().requireFIPNUM(); } + + const Grid& grid() + { return ebosSimulator_.gridManager().grid(); } + + const Schedule& schedule() const + { return ebosSimulator_.gridManager().schedule(); } + + const SummaryConfig& summaryConfig() const + { return ebosSimulator_.gridManager().summaryConfig(); } + + const EclipseState& eclState() const + { return ebosSimulator_.gridManager().eclState(); } + + bool isRestart() const { + const auto& initconfig = eclState().getInitConfig(); + 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 + 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 parallelOutput_; + const bool restart_double_si_; + std::unique_ptr< ThreadHandle > asyncOutput_; + }; + + + + + +} +#endif diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index bc4205429..cdfda731b 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -539,6 +539,8 @@ namespace Opm schedule(), summaryConfig())); eclIO_->writeInitial(computeLegacySimProps_(), int_vectors, nnc_); + Problem& problem = ebosProblem(); + problem.setEclIO(std::move(eclIO_)); } } @@ -550,13 +552,10 @@ namespace Opm // create output writer after grid is distributed, otherwise the parallel output // won't work correctly since we need to create a mapping from the distributed to // the global view - output_writer_.reset(new OutputWriter(grid(), - param_, - eclState(), - schedule(), - summaryConfig(), - std::move(eclIO_), - Opm::phaseUsageFromDeck(deck())) ); + + output_writer_.reset(new OutputWriter(*ebosSimulator_, + param_)); + } // Run the simulator. diff --git a/opm/autodiff/ParallelDebugOutput.hpp b/opm/autodiff/ParallelDebugOutput.hpp index d1d332920..61d9a336b 100644 --- a/opm/autodiff/ParallelDebugOutput.hpp +++ b/opm/autodiff/ParallelDebugOutput.hpp @@ -30,6 +30,7 @@ #include #include +#include #include #if HAVE_OPM_GRID diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 9964bea19..ad4094ac2 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -22,7 +22,7 @@ #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED -#include +#include #include #include #include @@ -61,7 +61,7 @@ public: typedef WellStateFullyImplicitBlackoil WellState; typedef BlackoilState ReservoirState; - typedef BlackoilOutputWriter OutputWriter; + typedef BlackoilOutputEbos OutputWriter; typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; @@ -140,20 +140,11 @@ public: failureReport_ = SimulatorReport(); if (output_writer_.isRestart()) { - // This is a restart, populate WellState and ReservoirState state objects from restart file + // This is a restart, populate WellState ReservoirState stateInit(Opm::UgGridHelpers::numCells(grid()), Opm::UgGridHelpers::numFaces(grid()), phaseUsage_.num_phases); output_writer_.initFromRestartFile(phaseUsage_, grid(), stateInit, prev_well_state, extra); - initHydroCarbonState(stateInit, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_); - initHysteresisParams(stateInit); - // communicate the restart solution to ebos - convertInput(/*iterationIdx=*/0, stateInit, ebosSimulator_ ); - ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); - // Sync the overlap region of the inital solution. It was generated - // from the ReservoirState which has wrong values in the ghost region - // for some models (SPE9, Norne, Model 2) - ebosSimulator_.model().syncOverlap(); } // Create timers and file for writing timing info. @@ -555,147 +546,6 @@ protected: const Schedule& schedule() const { return ebosSimulator_.gridManager().schedule(); } - void initHysteresisParams(ReservoirState& state) { - const int num_cells = Opm::UgGridHelpers::numCells(grid()); - - typedef std::vector VectorType; - - const VectorType& somax = state.getCellData( "SOMAX" ); - - for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) { - ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx); - } - - if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) { - auto matLawManager = ebosSimulator_.problem().materialLawManager(); - - VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" ); - VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" ); - - VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" ); - VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" ); - - for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) { - matLawManager->setOilWaterHysteresisParams( - pcSwMdc_ow[cellIdx], - krnSwMdc_ow[cellIdx], - cellIdx); - matLawManager->setGasOilHysteresisParams( - pcSwMdc_go[cellIdx], - krnSwMdc_go[cellIdx], - cellIdx); - } - } - } - - - // Used to convert initial Reservoirstate to primary variables in the SolutionVector - void convertInput( const int iterationIdx, - const ReservoirState& reservoirState, - Simulator& simulator ) const - { - SolutionVector& solution = simulator.model().solution( 0 /* timeIdx */ ); - const Opm::PhaseUsage pu = phaseUsage_; - - const std::vector active = detail::activePhases(pu); - bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); - bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); - - const int numCells = reservoirState.numCells(); - const int numPhases = phaseUsage_.num_phases; - const auto& oilPressure = reservoirState.pressure(); - const auto& saturations = reservoirState.saturation(); - const auto& rs = reservoirState.gasoilratio(); - const auto& rv = reservoirState.rv(); - for( int cellIdx = 0; cellIdx gas only with vaporized oil in the gas) is - // relatively expensive as it requires to compute the capillary - // pressure in order to get the gas phase pressure. (the reason why - // ebos uses the gas pressure here is that it makes the common case - // of the primary variable switching code fast because to determine - // whether the oil phase appears one needs to compute the Rv value - // for the saturated gas phase and if this is not available as a - // primary variable, it needs to be computed.) luckily for here, the - // gas-only case is not too common, so the performance impact of this - // is limited. - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - SatOnlyFluidState fluidState; - if ( active[Water] ) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); - } - else { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - } - fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); - - double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; - const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]); - - cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx]; - cellPv[BlackoilIndices::pressureSwitchIdx] = pg; - cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv ); - } - else - { - assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil); - cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]]; - cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ]; - cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg ); - } - } else { - // for oil-water case oil pressure should be used as primary variable - cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx]; - } - } - - // store the solution at the beginning of the time step - if( iterationIdx == 0 ) - { - simulator.model().solution( 1 /* timeIdx */ ) = solution; - } - } - // Data. Simulator& ebosSimulator_;