Merge remote-tracking branch 'upstream/master' into master-refactor-for-cpgrid-support

Conflicts:
	examples/sim_fibo_ad.cpp
	opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

examples/sim_2p_comp_ad.cpp

@@ -48,6 +48,8 @@
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 
 #include <boost/scoped_ptr.hpp>
 #include <boost/filesystem.hpp>
@@ -94,6 +96,7 @@ try
     boost::scoped_ptr<GridManager> grid;
     boost::scoped_ptr<BlackoilPropertiesInterface> props;
     boost::scoped_ptr<RockCompressibility> rock_comp;
+    EclipseStateConstPtr eclipseState;
     BlackoilState state;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -103,6 +106,8 @@ try
         deck.reset(new EclipseGridParser(deck_filename));
         Opm::ParserPtr newParser(new Opm::Parser() );
         Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename );
+        eclipseState.reset( new EclipseState(newParserDeck ));
+
 
         // Grid init
         grid.reset(new GridManager(newParserDeck));
@@ -255,7 +260,7 @@ try
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
             // Create new wells, well_state
-            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.

examples/sim_2p_incomp_ad.cpp

@@ -46,6 +46,9 @@
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
 
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 #include <boost/scoped_ptr.hpp>
 #include <boost/filesystem.hpp>
 
@@ -100,12 +103,16 @@ try
     boost::scoped_ptr<GridManager> grid;
     boost::scoped_ptr<IncompPropertiesInterface> props;
     boost::scoped_ptr<RockCompressibility> rock_comp;
+    EclipseStateConstPtr eclipseState;
     TwophaseState state;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     if (use_deck) {
         std::string deck_filename = param.get<std::string>("deck_filename");
+        ParserPtr parser(new Opm::Parser());
+        eclipseState.reset( new EclipseState(parser->parseFile(deck_filename)));
+
         deck.reset(new EclipseGridParser(deck_filename));
         // Grid init
         grid.reset(new GridManager(*deck));
@@ -262,7 +269,7 @@ try
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
             // Create new wells, well_state
-            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.

examples/sim_fibo_ad.cpp

@@ -236,22 +236,17 @@ try
             outputWriter.writeInit(simtimer, state, well_state.basicWellState());
 
         // Create and run simulator.
-    SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
+        SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
-                                                               *grid->c_grid(),
+                                                 *grid->c_grid(),
-                                                               *new_props,
+                                                 *new_props,
-                                                               rock_comp->isActive() ? rock_comp.get() : 0,
+                                                 rock_comp->isActive() ? rock_comp.get() : 0,
-                                                               wells,
+                                                 wells,
-                                                               linsolver,
+                                                 linsolver,
-                                                               grav,
+                                                 grav);
-                                                               outputWriter);
         SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
 
         outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
         fullReport += episodeReport;
-
-        if (output) {
-            episodeReport.reportParam(outStream);
-        }
     }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
@@ -318,8 +313,7 @@ try
                                                  rock_comp->isActive() ? rock_comp.get() : 0,
                                                  wells,
                                                  linsolver,
-                                                 grav,
+                                                 grav);
-                                                 outputWriter);
         outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
 
         if (epoch == 0) {

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -33,6 +33,7 @@
 #include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/Exceptions.hpp>
 #include <opm/core/well_controls.h>
 
 #ifdef HAVE_DUNE_CORNERPOINT
@@ -351,91 +352,19 @@ namespace {
     FullyImplicitBlackoilSolver<T>::constantState(const BlackoilState& x,
                                                   const WellStateFullyImplicitBlackoil&     xw)
     {
-        using namespace Opm::AutoDiffGrid;
+        auto state = variableState(x, xw);
-        const int nc = numCells(grid_);
-        const int np = x.numPhases();
 
-        // The block pattern assumes the following primary variables:
+        // HACK: throw away the derivatives. this may not be the most
-        //    pressure
+        // performant way to do things, but it will make the state
-        //    water saturation (if water present)
+        // automatically consistent with variableState() (and doing
-        //    gas saturation, Rv (vapor oil/gas ratio) or Rs (solution gas/oil ratio) depending on hydrocarbon state
+        // things automatically is all the rage in this module ;)
-        //      Gas only (undersaturated gas): Rv
+        state.pressure = ADB::constant(state.pressure.value());
-        //      Gas and oil: Sg
+        state.rs = ADB::constant(state.rs.value());
-        //      Oil only (undersaturated oil): Rs
+        state.rv = ADB::constant(state.rv.value());
-        //    well rates per active phase and well
+        for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx)
-        //    well bottom-hole pressure
+            state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
-        // Note that oil is assumed to always be present, but is never
+        state.qs = ADB::constant(state.qs.value());
-        // a primary variable.
+        state.bhp = ADB::constant(state.bhp.value());
-        assert(active_[ Oil ]);
-        std::vector<int> bpat(np, nc);
-        bpat.push_back(xw.bhp().size() * np);
-        bpat.push_back(xw.bhp().size());
-
-        SolutionState state(np);
-
-        // Pressure.
-        assert (not x.pressure().empty());
-        const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
-        state.pressure = ADB::constant(p, bpat);
-
-        // Saturation.
-        assert (not x.saturation().empty());
-        const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
-        const Opm::PhaseUsage pu = fluid_.phaseUsage();
-        {
-            V so = V::Ones(nc, 1);
-            if (active_[ Water ]) {
-                const int pos = pu.phase_pos[ Water ];
-                const V   sw  = s.col(pos);
-                so -= sw;
-
-                state.saturation[pos] = ADB::constant(sw, bpat);
-            }
-            if (active_[ Gas ]) {
-                const int pos = pu.phase_pos[ Gas ];
-                const V   sg  = s.col(pos);
-                so -= sg;
-
-                state.saturation[pos] = ADB::constant(sg, bpat);
-            }
-            if (active_[ Oil ]) {
-                const int pos = pu.phase_pos[ Oil ];
-                state.saturation[pos] = ADB::constant(so, bpat);
-            }
-        }
-
-        // Solution Gas-oil ratio (rs).
-        if (active_[ Oil ] && active_[ Gas ]) {
-            const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
-            state.rs = ADB::constant(rs, bpat);
-        } else {
-            const V Rs = V::Zero(nc, 1);
-            state.rs = ADB::constant(Rs, bpat);
-        }
-
-        // Vapor Oil-gas ratio (rv).
-        if (active_[ Oil ] && active_[ Gas ]) {
-            const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
-            state.rv = ADB::constant(rv, bpat);
-        } else {
-            const V rv = V::Zero(nc, 1);
-            state.rv = ADB::constant(rv, bpat);
-        }
-
-        // Well rates.
-        assert (not xw.wellRates().empty());
-        // Need to reshuffle well rates, from ordered by wells, then phase,
-        // to ordered by phase, then wells.
-        const int nw = wells_.number_of_wells;
-        // The transpose() below switches the ordering.
-        const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
-        const V qs = Eigen::Map<const V>(wrates.data(), nw*np);
-        state.qs = ADB::constant(qs, bpat);
-
-        // Well bottom-hole pressure.
-        assert (not xw.bhp().empty());
-        const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
-        state.bhp = ADB::constant(bhp, bpat);
 
         return state;
     }
@@ -1335,18 +1264,22 @@ void resetAnyComm(boost::any& anyComm, const Dune::CpGrid& grid)
     double
     FullyImplicitBlackoilSolver<T>::residualNorm() const
     {
-        double r = 0;
+        double globalNorm = 0;
-        for (std::vector<ADB>::const_iterator
+        std::vector<ADB>::const_iterator quantityIt = residual_.mass_balance.begin();
-                 b = residual_.mass_balance.begin(),
+        const std::vector<ADB>::const_iterator endQuantityIt = residual_.mass_balance.end();
-                 e = residual_.mass_balance.end();
+        for (; quantityIt != endQuantityIt; ++quantityIt)
-             b != e; ++b)
         {
-            r = std::max(r, (*b).value().matrix().norm());
+            const double quantityResid = (*quantityIt).value().matrix().norm();
+            if (!std::isfinite(quantityResid)) {
+                OPM_THROW(Opm::NumericalProblem,
+                          "Encountered a non-finite residual");
+            }
+            globalNorm = std::max(globalNorm, quantityResid);
         }
-        r = std::max(r, residual_.well_flux_eq.value().matrix().norm());
+        globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm());
-        r = std::max(r, residual_.well_eq.value().matrix().norm());
+        globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm());
 
-        return r;
+        return globalNorm;
     }
 
 

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -31,7 +31,6 @@ namespace Opm
 {
     namespace parameter { class ParameterGroup; }
     class BlackoilPropsAdInterface;
-    class EclipseWriter;
     class RockCompressibility;
     class WellsManager;
     class LinearSolverInterface;
@@ -75,8 +74,7 @@ namespace Opm
                                        const RockCompressibility* rock_comp_props,
                                        WellsManager& wells_manager,
                                        LinearSolverInterface& linsolver,
-                                       const double* gravity,
+                                       const double* gravity);
-                                       EclipseWriter &writer);
 
         /// Run the simulation.
         /// This will run succesive timesteps until timer.done() is true. It will

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -33,7 +33,6 @@
 #include <opm/core/simulator/SimulatorReport.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
-#include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
@@ -67,8 +66,7 @@ namespace Opm
              const RockCompressibility* rock_comp_props,
              WellsManager& wells_manager,
              LinearSolverInterface& linsolver,
-             const double* gravity,
+             const double* gravity);
-             EclipseWriter &writer);
 
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
@@ -97,7 +95,6 @@ namespace Opm
         FullyImplicitBlackoilSolver<Grid> solver_;
         // Misc. data
         std::vector<int> allcells_;
-        EclipseWriter &eclipseWriter_;
     };
 
 
@@ -110,11 +107,10 @@ namespace Opm
                                                                    const RockCompressibility* rock_comp_props,
                                                                    WellsManager& wells_manager,
                                                                    LinearSolverInterface& linsolver,
-                                                                   const double* gravity,
+                                                                   const double* gravity)
-                                                                   EclipseWriter &eclipseWriter)
 
     {
-        pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, eclipseWriter));
+        pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity));
     }
 
 
@@ -196,8 +192,7 @@ namespace Opm
                                                const RockCompressibility* rock_comp_props,
                                                WellsManager& wells_manager,
                                                LinearSolverInterface& linsolver,
-                                               const double* gravity,
+                                               const double* gravity)
-                                               EclipseWriter &eclipseWriter)
         : grid_(grid),
           props_(props),
           rock_comp_props_(rock_comp_props),
@@ -205,9 +200,7 @@ namespace Opm
           wells_(wells_manager.c_wells()),
           gravity_(gravity),
           geo_(grid_, props_, gravity_),
-          solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver),
+          solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
-          eclipseWriter_(eclipseWriter)
-
           /*                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
                                param.getDefault("nl_pressure_change_tolerance", 1.0),
                                param.getDefault("nl_pressure_maxiter", 10),
@@ -351,11 +344,6 @@ namespace Opm
 
             // advance to next timestep before reporting at this location
             ++timer;
-
-            // write an output file for later inspection
-            if (output_) {
-                eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
-            }
         }
 
         total_timer.stop();