Merge pull request #542 from andlaus/implement_MULTFLT_in_SCHEDULE_section

implement abililty to change the geology during the simulation run

examples/flow.cpp

@@ -434,7 +434,7 @@ try
                       << std::flush;
         }
 
-        SimulatorReport fullReport = simulator.run(simtimer, state);
+        SimulatorReport fullReport = simulator.run(eclipseState, simtimer, state);
 
         if( output_cout )
         {

examples/flow_solvent.cpp

@@ -407,7 +407,7 @@ try
                       << std::flush;
         }
 
-        SimulatorReport fullReport = simulator.run(simtimer, state);
+        SimulatorReport fullReport = simulator.run(eclipseState, simtimer, state);
 
         if( output_cout )
         {

opm/autodiff/GeoProps.hpp

@@ -69,6 +69,17 @@ namespace Opm
             , trans_(Opm::AutoDiffGrid::numFaces(grid))
             , gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
             , z_(Opm::AutoDiffGrid::numCells(grid))
+            , use_local_perm_(use_local_perm)
+        {
+            update(grid, props, eclState, grav);
+        }
+
+        /// compute the all geological properties at a given report step
+        template <class Props, class Grid>
+        void update(const Grid&              grid,
+                    const Props&             props ,
+                    Opm::EclipseStateConstPtr eclState,
+                    const double*            grav)
         {
             int numCells = AutoDiffGrid::numCells(grid);
             int numFaces = AutoDiffGrid::numFaces(grid);
@@ -121,7 +132,7 @@ namespace Opm
             Vector htrans(AutoDiffGrid::numCellFaces(grid));
             Grid* ug = const_cast<Grid*>(& grid);
 
-            if (! use_local_perm) {
+            if (! use_local_perm_) {
                 tpfa_htrans_compute(ug, props.permeability(), htrans.data());
             }
             else {
@@ -205,9 +216,7 @@ namespace Opm
         Vector gpot_ ;
         Vector z_;
         double gravity_[3]; // Size 3 even if grid is 2-dim.
-
+        bool use_local_perm_;
-
-
 
     };
 

opm/autodiff/SimulatorBase.hpp

@@ -120,7 +120,7 @@ namespace Opm
         /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
         SimulatorBase(const parameter::ParameterGroup& param,
                       const Grid& grid,
-                      const DerivedGeology& geo,
+                      DerivedGeology& geo,
                       BlackoilPropsAdInterface& props,
                       const RockCompressibility* rock_comp_props,
                       NewtonIterationBlackoilInterface& linsolver,
@@ -134,11 +134,13 @@ namespace Opm
         /// Run the simulation.
         /// This will run succesive timesteps until timer.done() is true. It will
         /// modify the reservoir and well states.
+        /// \param[in] eclState        the object which represents an internalized ECL deck
         /// \param[in,out] timer       governs the requested reporting timesteps
         /// \param[in,out] state       state of reservoir: pressure, fluxes
         /// \param[in,out] well_state  state of wells: bhp, perforation rates
         /// \return                    simulation report, with timing data
-        SimulatorReport run(SimulatorTimer& timer,
+        SimulatorReport run(EclipseStateConstPtr eclState,
+                            SimulatorTimer& timer,
                             ReservoirState& state);
 
     protected:
@@ -176,7 +178,7 @@ namespace Opm
         const RockCompressibility* rock_comp_props_;
         const double* gravity_;
         // Solvers
-        const DerivedGeology& geo_;
+        DerivedGeology& geo_;
         NewtonIterationBlackoilInterface& solver_;
         // Misc. data
         std::vector<int> allcells_;

opm/autodiff/SimulatorBase_impl.hpp

@@ -27,7 +27,7 @@ namespace Opm
     template <class Implementation>
     SimulatorBase<Implementation>::SimulatorBase(const parameter::ParameterGroup& param,
                                                  const Grid& grid,
-                                                 const DerivedGeology& geo,
+                                                 DerivedGeology& geo,
                                                  BlackoilPropsAdInterface& props,
                                                  const RockCompressibility* rock_comp_props,
                                                  NewtonIterationBlackoilInterface& linsolver,
@@ -74,7 +74,8 @@ namespace Opm
     }
 
     template <class Implementation>
-    SimulatorReport SimulatorBase<Implementation>::run(SimulatorTimer& timer,
+    SimulatorReport SimulatorBase<Implementation>::run(EclipseStateConstPtr eclState,
+                                                       SimulatorTimer& timer,
                                                        ReservoirState& state)
     {
         WellState prev_well_state;
@@ -88,6 +89,9 @@ namespace Opm
         std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
         std::ofstream tstep_os(tstep_filename.c_str());
 
+        const auto& schedule = eclState->getSchedule();
+        const auto& events = schedule->getEvents();
+
         // adaptive time stepping
         std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
         if( param_.getDefault("timestep.adaptive", true ) )
@@ -164,6 +168,21 @@ namespace Opm
                 solver->step(timer.currentStepLength(), state, well_state);
             }
 
+            // update the derived geology (transmissibilities, pore volumes, etc) if the
+            // has geology changed for the next report step
+            const int nextTimeStepIdx = timer.currentStepNum() + 1;
+            if (nextTimeStepIdx < timer.numSteps()
+                && events.hasEvent(ScheduleEvents::GEO_MODIFIER, nextTimeStepIdx)) {
+                // bring the contents of the keywords to the current state of the SCHEDULE
+                // section
+                //
+                // TODO (?): handle the parallel case (maybe this works out of the box)
+                ScheduleConstPtr schedule = eclipse_state_->getSchedule();
+                DeckConstPtr miniDeck = schedule->getModifierDeck(nextTimeStepIdx);
+                eclipse_state_->applyModifierDeck(miniDeck);
+                geo_.update(grid_, props_, eclipse_state_, gravity_);
+            }
+
             // take time that was used to solve system for this reportStep
             solver_timer.stop();
 

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -50,7 +50,7 @@ public:
     // forward the constructor to the base class
     SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
                                    const typename Base::Grid& grid,
-                                   const DerivedGeology& geo,
+                                   DerivedGeology& geo,
                                    BlackoilPropsAdInterface& props,
                                    const RockCompressibility* rock_comp_props,
                                    NewtonIterationBlackoilInterface& linsolver,

opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp

@@ -106,7 +106,7 @@ namespace Opm
     public:
         SimulatorFullyImplicitBlackoilSolvent(const parameter::ParameterGroup& param,
                                               const GridT& grid,
-                                              const DerivedGeology& geo,
+                                              DerivedGeology& geo,
                                               BlackoilPropsAdInterface& props,
                                               const SolventPropsAdFromDeck& solvent_props,
                                               const RockCompressibility* rock_comp_props,

opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp

@@ -26,7 +26,7 @@ namespace Opm
     SimulatorFullyImplicitBlackoilSolvent<GridT>::
     SimulatorFullyImplicitBlackoilSolvent(const parameter::ParameterGroup& param,
                                           const GridT& grid,
-                                          const DerivedGeology& geo,
+                                          DerivedGeology& geo,
                                           BlackoilPropsAdInterface& props,
                                           const SolventPropsAdFromDeck& solvent_props,
                                           const RockCompressibility* rock_comp_props,