Mearge from upstream

examples/Makefile.am

@@ -4,7 +4,10 @@ AM_CPPFLAGS =					\
 $(BOOST_CPPFLAGS)
 
 # All targets link to the library
-LDADD = $(top_builddir)/libopmcore.la
+LDADD =						\
+$(top_builddir)/libopmcore.la			\
+$(BOOST_FILESYSTEM_LIB)				\
+$(BOOST_SYSTEM_LIB)
 
 # ----------------------------------------------------------------------
 # Declare products (i.e., the example programs).
@@ -37,4 +40,7 @@ if UMFPACK
 noinst_PROGRAMS += spu_2p
 
 spu_2p_SOURCES   = spu_2p.cpp
+spu_2p_LDADD     =				\
+$(LDADD)					\
+$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 endif

examples/sim_2p_incomp_reorder.cpp

@@ -79,43 +79,34 @@ main(int argc, char** argv)
 
     // If we have a "deck_filename", grid and props will be read from that.
     bool use_deck = param.has("deck_filename");
+    boost::scoped_ptr<Opm::EclipseGridParser> deck;
     boost::scoped_ptr<Opm::GridManager> grid;
     boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
-    boost::scoped_ptr<Opm::WellsManager> wells;
     boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
-    Opm::SimulatorTimer simtimer;
     Opm::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");
-        Opm::EclipseGridParser deck(deck_filename);
+        deck.reset(new Opm::EclipseGridParser(deck_filename));
         // Grid init
-        grid.reset(new Opm::GridManager(deck));
+        grid.reset(new Opm::GridManager(*deck));
         // Rock and fluid init
         const int* gc = grid->c_grid()->global_cell;
         std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
-        props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
+        props.reset(new Opm::IncompPropertiesFromDeck(*deck, global_cell));
-        // Wells init.
-        wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
         // check_well_controls = param.getDefault("check_well_controls", false);
         // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
-        // Timer init.
-        if (deck.hasField("TSTEP")) {
-            simtimer.init(deck);
-        } else {
-            simtimer.init(param);
-        }
         // Rock compressibility.
-        rock_comp.reset(new Opm::RockCompressibility(deck));
+        rock_comp.reset(new Opm::RockCompressibility(*deck));
         // Gravity.
-        gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
+        gravity[2] = deck->hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
         // Init state variables (saturation and pressure).
         if (param.has("init_saturation")) {
             initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
         } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+            initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
         }
     } else {
         // Grid init.
@@ -128,10 +119,6 @@ main(int argc, char** argv)
         grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
         // Rock and fluid init.
         props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Wells init.
-        wells.reset(new Opm::WellsManager());
-        // Timer init.
-        simtimer.init(param);
         // Rock compressibility.
         rock_comp.reset(new Opm::RockCompressibility(param));
         // Gravity.
@@ -165,9 +152,8 @@ main(int argc, char** argv)
 
     // Initialising src
     std::vector<double> src(num_cells, 0.0);
-    if (wells->c_wells()) {
+    if (use_deck) {
         // Do nothing, wells will be the driving force, not source terms.
-        // Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
     } else {
         const double default_injection = use_gravity ? 0.0 : 0.1;
         const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
@@ -190,30 +176,69 @@ main(int argc, char** argv)
     const double *grav = use_gravity ? &gravity[0] : 0;
 
 
-    Opm::SimulatorTwophase simulator(param,
-                                     *grid->c_grid(),
-                                     *props,
-                                     rock_comp->isActive() ? rock_comp.get() : 0,
-                                     wells->c_wells(),
-                                     src,
-                                     bcs.c_bcs(),
-                                     linsolver,
-                                     grav);
-
     // Warn if any parameters are unused.
-    if (param.anyUnused()) {
+    // if (param.anyUnused()) {
-        std::cout << "--------------------   Unused parameters:   --------------------\n";
+    //     std::cout << "--------------------   Unused parameters:   --------------------\n";
-        param.displayUsage();
+    //     param.displayUsage();
-        std::cout << "----------------------------------------------------------------" << std::endl;
+    //     std::cout << "----------------------------------------------------------------" << std::endl;
-    }
+    // }
 
     // Write parameters used for later reference.
     // if (output) {
     //     param.writeParam(output_dir + "/spu_2p.param");
     // }
 
-    WellState well_state;
-    well_state.init(wells->c_wells(), state);
 
-    simulator.run(simtimer, state, well_state);
+    if (!use_deck) {
+        // Simple simulation without a deck.
+        Opm::SimulatorTwophase simulator(param,
+                                         *grid->c_grid(),
+                                         *props,
+                                         rock_comp->isActive() ? rock_comp.get() : 0,
+                                         0, // wells
+                                         src,
+                                         bcs.c_bcs(),
+                                         linsolver,
+                                         grav);
+        Opm::SimulatorTimer simtimer;
+        simtimer.init(param);
+        WellState well_state;
+        well_state.init(0, state);
+        simulator.run(simtimer, state, well_state);
+    } else {
+        // With a deck, we may have more epochs etc.
+        WellState well_state;
+        int step = 0;
+        for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
+            deck->setCurrentEpoch(epoch);
+            Opm::WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            Opm::SimulatorTwophase simulator(param,
+                                             *grid->c_grid(),
+                                             *props,
+                                             rock_comp->isActive() ? rock_comp.get() : 0,
+                                             wells.c_wells(),
+                                             src,
+                                             bcs.c_bcs(),
+                                             linsolver,
+                                             grav);
+            // @@@ 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.
+            if (epoch == 0) {
+                well_state.init(wells.c_wells(), state);
+            }
+            Opm::SimulatorTimer simtimer;
+            if (deck->hasField("TSTEP")) {
+                simtimer.init(*deck);
+            } else {
+                if (epoch != 0) {
+                    THROW("No TSTEP in deck for epoch " << epoch);
+                }
+                simtimer.init(param);
+            }
+            simtimer.setCurrentStepNum(step);
+            simulator.run(simtimer, state, well_state);
+            step = simtimer.currentStepNum();
+        }
+    }
 }

examples/sim_wateroil.cpp

@@ -273,12 +273,9 @@ main(int argc, char** argv)
     }
     double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
 
-    // We need a separate reorder_sat, because the reorder
-    // code expects a scalar sw, not both sw and so.
-    std::vector<double> reorder_sat(num_cells);
-    std::vector<double> src(num_cells, 0.0);
 
     // Initialising src
+    std::vector<double> src(num_cells, 0.0);
     if (wells->c_wells()) {
         // Do nothing, wells will be the driving force, not source terms.
         // Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
@@ -426,15 +423,12 @@ main(int argc, char** argv)
             std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl;
         }
         for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
-            Opm::toWaterSat(state.saturation(), reorder_sat);
             reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0],
-                                &porevol[0], &reorder_src[0], stepsize, &reorder_sat[0]);
+                                &porevol[0], &reorder_src[0], stepsize, state.saturation());
-            Opm::toBothSat(reorder_sat, state.saturation());
             // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
             if (use_segregation_split) {
                 THROW("Segregation not implemented yet.");
-                // reorder_model.solveGravity(columns, &porevol[0], stepsize, reorder_sat);
+                // reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
-                Opm::toBothSat(reorder_sat, state.saturation());
             }
         }
 #endif // TRANSPORT_SOLVER_FIXED

opm/core/eclipse/EclipseGridParser.cpp

@@ -395,7 +395,6 @@ void EclipseGridParser::readImpl(istream& is)
                     // Set it to START date, or default if no START.
                     // This will only ever happen in the first epoch,
                     // upon first encountering a timestepping keyword.
-                    SpecialMap::const_iterator it = sm.find("START");
                     if (hasField("START")) {
                         start_date_ = getSTART().date;
                     } else {
@@ -446,12 +445,27 @@ void EclipseGridParser::readImpl(istream& is)
                     special_field_by_epoch_.push_back(SpecialMap());
                     ++current_epoch_;
                     ASSERT(int(special_field_by_epoch_.size()) == current_epoch_ + 1);
+                    // Add clones of all existing special fields to new map.
+                    SpecialMap& oldmap = special_field_by_epoch_[current_epoch_ - 1];
+                    SpecialMap& newmap = special_field_by_epoch_[current_epoch_];
+                    for (SpecialMap::iterator it = oldmap.begin(); it != oldmap.end(); ++it) {
+                        // if (it->first != "TSTEP") {
+                            newmap[it->first] = cloneSpecialField(it->first, it->second);
+                            //}
+                    }
+                    //ASSERT(newmap.count("TSTEP") == 0);
                 }
-                SpecialFieldPtr sb_ptr = createSpecialField(is, keyword);
+                // Check if the keyword already exists. If so, append. Otherwise, create new.
-                if (sb_ptr) {
+                SpecialMap::iterator it = special_field_by_epoch_[current_epoch_].find(keyword);
-                    special_field_by_epoch_[current_epoch_][keyword] = sb_ptr;
+                if (it != special_field_by_epoch_[current_epoch_].end()) {
+                    it->second->read(is);
                 } else {
-                    THROW("Could not create field " << keyword);
+                    SpecialFieldPtr sb_ptr = createSpecialField(is, keyword);
+                    if (sb_ptr) {
+                        special_field_by_epoch_[current_epoch_][keyword] = sb_ptr;
+                    } else {
+                        THROW("Could not create field " << keyword);
+                    }
                 }
                 break;
             }
@@ -697,6 +711,19 @@ EclipseGridParser::createSpecialField(std::istream& is,
     return spec_ptr;
 }
 
+//---------------------------------------------------------------------------
+std::tr1::shared_ptr<SpecialBase>
+EclipseGridParser::cloneSpecialField(const std::string& fieldname,
+                                     const std::tr1::shared_ptr<SpecialBase> original)
+//---------------------------------------------------------------------------
+{
+    string ukey = upcase(fieldname);
+    std::tr1::shared_ptr<SpecialBase> spec_ptr
+        = Factory<SpecialBase>::cloneObject(fieldname, original);
+    return spec_ptr;
+}
+
+
 //---------------------------------------------------------------------------
 void EclipseGridParser::setIntegerField(const std::string& keyword,
                                         const std::vector<int>& field)

opm/core/eclipse/EclipseGridParser.hpp

@@ -193,6 +193,8 @@ public:
 
 private:
     SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname);
+    SpecialFieldPtr cloneSpecialField(const std::string& fieldname,
+                                      const std::tr1::shared_ptr<SpecialBase> original);
     void readImpl(std::istream& is);
 
 

opm/core/pressure/IncompTpfa.hpp

@@ -57,7 +57,7 @@ namespace Opm
 		   const double* permeability,
 		   const double* gravity,
                    LinearSolverInterface& linsolver,
-                   const struct Wells* wells = 0);
+                   const Wells* wells);
 
 	/// Destructor.
 	~IncompTpfa();

opm/core/simulator/SimulatorTwophase.cpp

@@ -214,7 +214,7 @@ namespace Opm
           bcs_(bcs),
           linsolver_(linsolver),
           gravity_(gravity),
-          psolver_(grid, props.permeability(), gravity, linsolver),
+          psolver_(grid, props.permeability(), gravity, linsolver, wells),
           tsolver_(grid, props, 1e-9, 30)
     {
         // For output.

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -49,7 +49,7 @@ namespace Opm
 	  darcyflux_(0),
 	  source_(0),
 	  dt_(0.0),
-	  saturation_(0),
+	  saturation_(grid.number_of_cells, -1.0),
 	  fractionalflow_(grid.number_of_cells, -1.0),
 	  mob_(2*grid.number_of_cells, -1.0),
           ia_upw_(grid.number_of_cells + 1, -1),
@@ -79,14 +79,15 @@ namespace Opm
                                                    const double* porevolume,
                                                    const double* source,
                                                    const double dt,
-                                                   double* saturation)
+                                                   std::vector<double>& saturation)
     {
 	darcyflux_ = darcyflux;
         surfacevol0_ = surfacevol0;
         porevolume_ = porevolume;
 	source_ = source;
 	dt_ = dt;
-	saturation_ = saturation;
+        toWaterSat(saturation, saturation_);
+
         props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
         props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
 
@@ -103,6 +104,7 @@ namespace Opm
 			       &seq[0], &comp[0], &ncomp,
 			       &ia_downw_[0], &ja_downw_[0]);
 	reorderAndTransport(grid_, darcyflux);
+        toBothSat(saturation_, saturation);
     }
 
     // Residual function r(s) for a single-cell implicit Euler transport
@@ -489,7 +491,7 @@ namespace Opm
         // Set up other variables.
         porevolume_ = porevolume;
         dt_ = dt;
-        saturation_ = &saturation[0];
+        toWaterSat(saturation, saturation_);
 
         // Solve on all columns.
         int num_iters = 0;
@@ -499,6 +501,7 @@ namespace Opm
         }
         std::cout << "Gauss-Seidel column solver average iterations: "
                   << double(num_iters)/double(columns.size()) << std::endl;
+        toBothSat(saturation_, saturation);
     }
 
 } // namespace Opm

opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp

@@ -33,33 +33,56 @@ namespace Opm
     class TransportModelCompressibleTwophase : public TransportModelInterface
     {
     public:
+	/// Construct solver.
+	/// \param[in] grid      A 2d or 3d grid.
+	/// \param[in] props     Rock and fluid properties.
+	/// \param[in] tol       Tolerance used in the solver.
+	/// \param[in] maxit     Maximum number of non-linear iterations used.
 	TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
                                            const Opm::BlackoilPropertiesInterface& props,
                                            const double tol,
                                            const int maxit);
 
+	/// Solve for saturation at next timestep.
+	/// \param[in] darcyflux         Array of signed face fluxes.
+	/// \param[in] pressure          Array of cell pressures
+	/// \param[in] surfacevol0       Array of surface volumes at start of timestep
+	/// \param[in] porevolume        Array of pore volumes.
+	/// \param[in] source            Transport source term.
+	/// \param[in] dt                Time step.
+	/// \param[in, out] saturation   Phase saturations.
 	void solve(const double* darcyflux,
                    const double* pressure,
                    const double* surfacevol0,
                    const double* porevolume,
 		   const double* source,
 		   const double dt,
-		   double* saturation);
+		   std::vector<double>& saturation);
-
-	virtual void solveSingleCell(const int cell);
-	virtual void solveMultiCell(const int num_cells, const int* cells);
 
+        /// Initialise quantities needed by gravity solver.
+        /// \param[in] grav    gravity vector
         void initGravity(const double* grav);
-        void solveSingleCellGravity(const std::vector<int>& cells,
+
-                                    const int pos,
+        /// Solve for gravity segregation.
-                                    const double* gravflux);
+        /// This uses a column-wise nonlinear Gauss-Seidel approach.
-        int solveGravityColumn(const std::vector<int>& cells);
+        /// It assumes that the input columns contain cells in a single
+        /// vertical stack, that do not interact with other columns (for
+        /// gravity segregation.
+        /// \TODO: Implement this.
         void solveGravity(const std::vector<std::vector<int> >& columns,
                           const double* pressure,
                           const double* porevolume,
                           const double dt,
                           std::vector<double>& saturation);
 
+    private:
+	virtual void solveSingleCell(const int cell);
+	virtual void solveMultiCell(const int num_cells, const int* cells);
+        void solveSingleCellGravity(const std::vector<int>& cells,
+                                    const int pos,
+                                    const double* gravflux);
+        int solveGravityColumn(const std::vector<int>& cells);
+
     private:
 	const UnstructuredGrid& grid_;
 	const BlackoilPropertiesInterface& props_;
@@ -76,7 +99,7 @@ namespace Opm
 	const double* porevolume_;  // one volume per cell
 	const double* source_;      // one source per cell
 	double dt_;
-	double* saturation_;        // one per cell
+        std::vector<double> saturation_;        // P (= num. phases) per cell
 	std::vector<double> fractionalflow_;  // = m[0]/(m[0] + m[1]) per cell
         // For gravity segregation.
         std::vector<double> gravflux_;

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -57,8 +57,19 @@ namespace Opm
 		   const double dt,
 		   std::vector<double>& saturation);
 
+        /// Initialise quantities needed by gravity solver.
+        /// \param[in] grav    gravity vector
         void initGravity(const double* grav);
 
+        /// Solve for gravity segregation.
+        /// This uses a column-wise nonlinear Gauss-Seidel approach.
+        /// It assumes that the input columns contain cells in a single
+        /// vertical stack, that do not interact with other columns (for
+        /// gravity segregation.
+	/// \param[in] columns           Vector of cell-columns.
+	/// \param[in] porevolume        Array of pore volumes.
+	/// \param[in] dt                Time step.
+	/// \param[in, out] saturation   Phase saturations.
         void solveGravity(const std::vector<std::vector<int> >& columns,
                           const double* porevolume,
                           const double dt,

opm/core/utility/Factory.hpp

@@ -65,6 +65,18 @@ namespace Opm
 	    return instance().doCreateObject(type);
 	}
 
+	/// Clones an new object of the class associated with the given type string,
+        /// and returns a pointer to it.
+	/// \param type the type string of the class that the user wants to have
+	///             constructed.
+	/// \param original (smart) pointer to object to be cloned.
+        /// \return (smart) pointer to the created object.
+	static ProductPtr cloneObject(const std::string& type,
+                                      const ProductPtr original)
+	{
+	    return instance().doCloneObject(type, original);
+	}
+
 	/// Add a creator to the Factory.
         /// After the call, the user may obtain new objects of the Derived type by
         /// calling createObject() with the given type string as an argument.
@@ -97,6 +109,7 @@ namespace Opm
         {
         public:
             virtual ProductPtr create() = 0;
+            virtual ProductPtr clone(ProductPtr original) = 0;
             virtual ~Creator() {}
         };
 
@@ -109,6 +122,11 @@ namespace Opm
             {
                 return ProductPtr(new Derived);
             }
+            virtual ProductPtr clone(const ProductPtr original)
+            {
+                const Derived& orig = dynamic_cast<const Derived&>(*original);
+                return ProductPtr(new Derived(orig));
+            }
         };
 
         typedef std::tr1::shared_ptr<Creator> CreatorPtr;
@@ -128,6 +146,19 @@ namespace Opm
 	    return it->second->create();
 	}
 
+        // Actually clones the product object.
+	ProductPtr doCloneObject(const std::string& type,
+                                 const ProductPtr original)
+	{
+	    typename CreatorMap::iterator it;
+	    it = string_to_creator_.find(type);
+	    if (it == string_to_creator_.end()) {
+		THROW("Creator type " << type
+		      << " is not registered in the factory.");
+	    }
+	    return it->second->clone(original);
+	}
+
         // Actually adds the creator.
         template <class Derived>
         void doAddCreator(const std::string& type)

opm/core/utility/SimulatorTimer.cpp

@@ -66,6 +66,13 @@ namespace Opm
 	return current_step_;
     }
 
+    /// Set current step number.
+    void SimulatorTimer::setCurrentStepNum(int step)
+    {
+        current_step_ = step;
+    }
+
+
     /// Current step length.
     double SimulatorTimer::currentStepLength() const
     {

opm/core/utility/SimulatorTimer.hpp

@@ -50,6 +50,9 @@ namespace Opm
 	/// Current step number.
 	int currentStepNum() const;
 
+        /// Set current step number.
+        void setCurrentStepNum(int step);
+
 	/// Current step length.
 	/// Note: if done(), it is an error to call currentStepLength().
 	double currentStepLength() const;

tutorials/tutorial2.cpp

@@ -108,12 +108,15 @@ int main()
     Opm::LinearSolverUmfpack linsolver;
     /// \endcode
     /// \page tutorial2
+
     /// We set up a pressure solver for the incompressible problem,
-    /// using the two-point flux approximation discretization.
+    /// using the two-point flux approximation discretization.  The
-    /// The third argument which corresponds to gravity is set to
+    /// third argument which corresponds to gravity is set to a null
-    /// zero (no gravity).
+    /// pointer (no gravity). The final argument would be a pointer to
+    /// a Wells data structure, again we use a null pointer to
+    /// indicate that we have no wells.
     /// \code
-    Opm::IncompTpfa psolver(*grid.c_grid(), &permeability[0], 0, linsolver);
+    Opm::IncompTpfa psolver(*grid.c_grid(), &permeability[0], 0, linsolver, 0);
     /// \endcode
     /// \page tutorial2
     /// We define the source term.

tutorials/tutorial3.cpp

@@ -164,10 +164,11 @@ int main ()
     /// \page tutorial3
     /// \details We may now set up the pressure solver. At this point,
     /// unchanging parameters such as transmissibility are computed
-    /// and stored internally by the IncompTpfa class.
+    /// and stored internally by the IncompTpfa class. The final (null pointer)
+    /// constructor argument is for wells, which are now used in this tutorial.
     /// \code
     LinearSolverUmfpack linsolver;
-    IncompTpfa psolver(grid, props.permeability(), grav, linsolver);
+    IncompTpfa psolver(grid, props.permeability(), grav, linsolver, 0);
     /// \endcode
 
     /// \page tutorial3