diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 82675b09d..d660e3cf7 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -77,11 +77,6 @@ SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true); // default in flow is to formulate the equations in surface volumes SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true); SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual, false); - - -// SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy -// code for fluid and satfunc handling gets fully retired. -SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false); }} namespace Opm { @@ -1133,7 +1128,7 @@ namespace Opm { return regionValues; } - SimulationDataContainer getSimulatorData ( const SimulationDataContainer& localState) const + SimulationDataContainer getSimulatorData ( const SimulationDataContainer& /*localState*/) const { typedef std::vector VectorType; @@ -1341,43 +1336,37 @@ namespace Opm { // 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 passed by the localState. + // 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 ) { - Rs[cellIdx] = localState.getCellData( BlackoilState::GASOILRATIO )[cellIdx]; - Rv[cellIdx] = localState.getCellData( BlackoilState::RV)[cellIdx]; + typedef Opm::CompositionalFluidState ScalarFluidState; + const ScalarFluidState& fs_updated = ebosSimulator().problem().initialFluidState(cellIdx); - // copy the fluidstate and set the new rs and rv values - auto fs_updated = fs; - auto rs_eval = fs_updated.Rs(); - rs_eval.setValue( Rs[cellIdx] ); - fs_updated.setRs(rs_eval); - auto rv_eval = fs_updated.Rv(); - rv_eval.setValue( Rv[cellIdx] ); - fs_updated.setRv(rv_eval); + // use initial rs and rv values + Rv[cellIdx] = Opm::BlackOil::getRv_(fs_updated, intQuants.pvtRegionIndex()); + Rs[cellIdx] = Opm::BlackOil::getRs_(fs_updated, intQuants.pvtRegionIndex()); //re-compute the volume factors, viscosities and densities. rhoOil[cellIdx] = FluidSystem::density(fs_updated, FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()).value(); + intQuants.pvtRegionIndex()); rhoGas[cellIdx] = FluidSystem::density(fs_updated, FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()).value(); + intQuants.pvtRegionIndex()); bOil[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated, FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()).value(); + intQuants.pvtRegionIndex()); bGas[cellIdx] = FluidSystem::inverseFormationVolumeFactor(fs_updated, FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()).value(); - + intQuants.pvtRegionIndex()); muOil[cellIdx] = FluidSystem::viscosity(fs_updated, FluidSystem::oilPhaseIdx, - intQuants.pvtRegionIndex()).value(); + intQuants.pvtRegionIndex()); muGas[cellIdx] = FluidSystem::viscosity(fs_updated, FluidSystem::gasPhaseIdx, - intQuants.pvtRegionIndex()).value(); + intQuants.pvtRegionIndex()); } } diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 4158dc2bd..bc4205429 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -42,10 +42,6 @@ #include -#include -#include -#include - #include #include #include @@ -114,7 +110,6 @@ namespace Opm printPRTHeader(); extractMessages(); runDiagnostics(); - setupState(); writeInit(); setupOutputWriter(); setupLinearSolver(); @@ -481,118 +476,7 @@ namespace Opm const SummaryConfig& summaryConfig() const { return ebosSimulator_->gridManager().summaryConfig(); } - - // Initialise the reservoir state. Updated fluid props for SWATINIT. - // Writes to: - // state_ - // threshold_pressures_ - // fluidprops_ (if SWATINIT is used) - void setupState() - { - const PhaseUsage pu = Opm::phaseUsageFromDeck(deck()); - const Grid& grid = this->grid(); - - // Need old-style fluid object for init purposes (only). - BlackoilPropertiesFromDeck props(deck(), - eclState(), - materialLawManager(), - grid.size(/*codim=*/0), - grid.globalCell().data(), - grid.logicalCartesianSize().data(), - param_); - - - // Init state variables (saturation and pressure). - if (param_.has("init_saturation")) { - state_.reset(new ReservoirState(grid.size(/*codim=*/0), - grid.numFaces(), - props.numPhases())); - - initStateBasic(grid.size(/*codim=*/0), - grid.globalCell().data(), - grid.logicalCartesianSize().data(), - grid.numFaces(), - Opm::UgGridHelpers::faceCells(grid), - Opm::UgGridHelpers::beginFaceCentroids(grid), - Opm::UgGridHelpers::beginCellCentroids(grid), - Grid::dimension, - props, param_, gravity(), *state_); - - initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_); - - enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; - if (pu.phase_used[Oil] && pu.phase_used[Gas]) { - const int numPhases = props.numPhases(); - const int numCells = Opm::UgGridHelpers::numCells(grid); - - // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState. - auto& gor = state_->getCellData( BlackoilState::GASOILRATIO ); - const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL ); - for (int c = 0; c < numCells; ++c) { - // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field. - gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]]; - } - } - } else if (deck().hasKeyword("EQUIL")) { - // Which state class are we really using - what a f... mess? - state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::numFaces(grid), - props.numPhases())); - - typedef EQUIL::DeckDependent::InitialStateComputer ISC; - - ISC isc(*materialLawManager(), eclState(), grid, gravity()); - - const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const int oilpos = FluidSystem::oilPhaseIdx; - const int waterpos = FluidSystem::waterPhaseIdx; - const int ref_phase = oil ? oilpos : waterpos; - - state_->pressure() = isc.press()[ref_phase]; - convertSats(state_->saturation(), isc.saturation(), pu); - state_->gasoilratio() = isc.rs(); - state_->rv() = isc.rv(); - - //initStateEquil(grid, materialLawManager(), eclState(), gravity(), pu, *state_); - //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); - } else { - state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::numFaces(grid), - props.numPhases())); - initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::globalCell(grid), - Opm::UgGridHelpers::numFaces(grid), - Opm::UgGridHelpers::faceCells(grid), - Opm::UgGridHelpers::beginFaceCentroids(grid), - Opm::UgGridHelpers::beginCellCentroids(grid), - Opm::UgGridHelpers::dimensions(grid), - props, deck(), gravity(), *state_); - } - - initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL")); - - // Get initial polymer concentration from ebos - if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { - auto& cpolymer = state_->getCellData( state_->POLYMER ); - const int numCells = Opm::UgGridHelpers::numCells(grid); - for (int c = 0; c < numCells; ++c) { - cpolymer[c] = ebosProblem().polymerConcentration(c); - } - } - // Get initial solvent saturation from ebos - if (GET_PROP_VALUE(TypeTag, EnableSolvent)) { - auto& solvent = state_->getCellData( state_->SSOL ); - auto& sat = state_->saturation(); - const int np = props.numPhases(); - const int numCells = Opm::UgGridHelpers::numCells(grid); - for (int c = 0; c < numCells; ++c) { - solvent[c] = ebosProblem().solventSaturation(c); - sat[c * np + pu.phase_pos[Water]]; - } - } - - } - + // Extract messages from parser. // Writes to: // OpmLog singleton. @@ -695,7 +579,7 @@ namespace Opm OpmLog::info(msg); } - SimulatorReport successReport = simulator_->run(simtimer, *state_); + SimulatorReport successReport = simulator_->run(simtimer); SimulatorReport failureReport = simulator_->failureReport(); if (output_cout_) { @@ -830,6 +714,7 @@ namespace Opm typedef typename Grid::LeafGridView GridView; using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; // Get the owner rank number for each cell + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; using Handle = CellOwnerDataHandle; const Grid& globalGrid = this->globalGrid(); const auto& globalGridView = globalGrid.leafGridView(); @@ -1022,7 +907,6 @@ namespace Opm ParameterGroup param_; bool output_to_files_ = false; std::string output_dir_ = std::string("."); - std::unique_ptr state_; NNC nnc_; std::unique_ptr eclIO_; std::unique_ptr output_writer_; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 4928d6f6c..ce24dd942 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -130,36 +130,34 @@ public: /// \param[in,out] timer governs the requested reporting timesteps /// \param[in,out] state state of reservoir: pressure, fluxes /// \return simulation report, with timing data - SimulatorReport run(SimulatorTimer& timer, - ReservoirState& state) + SimulatorReport run(SimulatorTimer& timer) { + + ReservoirState dummy_state(0,0,0); + WellState prev_well_state; ExtraData extra; failureReport_ = SimulatorReport(); - // communicate the initial solution to ebos - if (timer.initialStep()) { - convertInput(/*iterationIdx=*/0, state, ebosSimulator_ ); - ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); - } - if (output_writer_.isRestart()) { // This is a restart, populate WellState and ReservoirState state objects from restart file - output_writer_.initFromRestartFile(phaseUsage_, grid(), state, prev_well_state, extra); - initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_); - initHysteresisParams(state); + 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, state, ebosSimulator_ ); + 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(); } - // 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. Opm::time::StopWatch solver_timer; Opm::time::StopWatch total_timer; @@ -227,7 +225,7 @@ public: // No per cell data is written for initial step, but will be // for subsequent steps, when we have started simulating - output_writer_.writeTimeStep( timer, state, well_model.wellState(), solver->model() ); + output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() ); report.output_write_time += perfTimer.stop(); } @@ -263,14 +261,14 @@ public: events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum()); - stepReport = adaptiveTimeStepping->step( timer, *solver, state, wellStateDummy, event, output_writer_, + stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, output_writer_.requireFIPNUM() ? &fipnum_ : nullptr ); report += stepReport; failureReport_ += adaptiveTimeStepping->failureReport(); } else { // solve for complete report step - stepReport = solver->step(timer, state, wellStateDummy); + stepReport = solver->step(timer, dummy_state, wellStateDummy); report += stepReport; failureReport_ += solver->failureReport(); @@ -298,13 +296,6 @@ public: // update timing. report.solver_time += solver_timer.secsSinceStart(); - // We don't need the reservoir state anymore. It is just passed around to avoid - // code duplication. Pass empty state instead. - if (timer.initialStep()) { - ReservoirState stateTrivial(0,0,0); - state = stateTrivial; - } - // Increment timer, remember well state. ++timer; @@ -325,7 +316,7 @@ public: Dune::Timer perfTimer; perfTimer.start(); const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0; - output_writer_.writeTimeStep( timer, 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(); }