Use the initial solution from Ebos

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<double> 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];
+                    typedef Opm::CompositionalFluidState<Scalar, FluidSystem> ScalarFluidState;
-                    Rv[cellIdx] = localState.getCellData( BlackoilState::RV)[cellIdx];
+                    const ScalarFluidState& fs_updated = ebosSimulator().problem().initialFluidState(cellIdx);
 
-                    // copy the fluidstate and set the new rs and rv values
+                    // use initial rs and rv values
-                    auto fs_updated = fs;
+                    Rv[cellIdx] = Opm::BlackOil::getRv_<FluidSystem, Scalar, ScalarFluidState>(fs_updated,  intQuants.pvtRegionIndex());
-                    auto rs_eval = fs_updated.Rs();
+                    Rs[cellIdx] = Opm::BlackOil::getRs_<FluidSystem, Scalar, ScalarFluidState>(fs_updated,  intQuants.pvtRegionIndex());
-                    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);
 
                     //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());
 
                 }
             }

opm/autodiff/FlowMainEbos.hpp

@@ -42,10 +42,6 @@
 
 #include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
 
-#include <opm/core/simulator/initStateEquil.hpp>
-#include <opm/core/simulator/initState.hpp>
-#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
-
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/OpmLog/EclipsePRTLog.hpp>
 #include <opm/common/OpmLog/LogUtil.hpp>
@@ -114,7 +110,6 @@ namespace Opm
                 printPRTHeader();
                 extractMessages();
                 runDiagnostics();
-                setupState();
                 writeInit();
                 setupOutputWriter();
                 setupLinearSolver();
@@ -482,117 +477,6 @@ 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<FluidSystem> 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<FluidSystem>(state_->saturation(), isc.saturation(), pu);
-                state_->gasoilratio() = isc.rs();
-                state_->rv() = isc.rv();
-
-                //initStateEquil<FluidSystem>(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<GridView, Dune::MCMGElementLayout>;
                 // Get the owner rank number for each cell
+                using ElementMapper =  Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout>;
                 using Handle = CellOwnerDataHandle<ElementMapper>;
                 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<ReservoirState> state_;
         NNC nnc_;
         std::unique_ptr<EclipseIO> eclIO_;
         std::unique_ptr<OutputWriter> output_writer_;

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,
+    SimulatorReport run(SimulatorTimer& timer)
-                        ReservoirState& state)
     {
+
+        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);
+            ReservoirState stateInit(Opm::UgGridHelpers::numCells(grid()),
-            initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
+                                     Opm::UgGridHelpers::numFaces(grid()),
-            initHysteresisParams(state);
+                                     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();
 
         }