Merge remote-tracking branch 'atgeirr/master'

CMakeLists_files.cmake

@@ -27,6 +27,7 @@
 # find opm -name '*.c*' -printf '\t%p\n' | sort
 list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/BlackoilPropsAd.cpp
+	opm/autodiff/BlackoilPropsAdInterface.cpp
 	opm/autodiff/SimulatorIncompTwophaseAdfi.cpp
 	opm/autodiff/TransportSolverTwophaseAd.cpp
 	)
@@ -65,6 +66,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/AutoDiffHelpers.hpp
 	opm/autodiff/AutoDiff.hpp
 	opm/autodiff/BlackoilPropsAd.hpp
+	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/ImpesTPFAAD.hpp
 	opm/autodiff/SimulatorIncompTwophaseAdfi.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp

examples/test_impestpfa_ad.cpp

@@ -119,20 +119,22 @@ main(int argc, char* argv[])
     typedef Opm::BlackoilPropertiesInterface    Geology;
     typedef DerivedGeology<Geology, ADB::V>     GeoProps;
     typedef Opm::BlackoilPropsAd    BOFluid;
-    typedef Opm::ImpesTPFAAD<BOFluid, GeoProps> PSolver;
+    typedef Opm::ImpesTPFAAD<GeoProps> PSolver;
 
     Wells* wells = create_wells(2, 2, 5);
     const double inj_frac[] = { 1.0, 0.0 };
     const double prod_frac[] = { 0.0, 0.0 };
-    const int inj_cells[] = { 0, 1, 2 };
+    const int num_inj = 3;
-    const int prod_cells[] = { 20, 21 };
+    const int inj_cells[num_inj] = { 0, 1, 2 };
-    const double WI[3] = { 1e-14 };
+    const int num_prod = 2;
-    bool ok = add_well(INJECTOR, 0.0, 1, inj_frac, inj_cells, WI, "Inj", wells);
+    const int prod_cells[num_prod] = { 20, 21 };
-    ok = ok && add_well(PRODUCER, 0.0, 1, prod_frac, prod_cells, WI, "Prod", wells);
+    const double WI[3] = { 1e-12, 1e-12, 1e-12 };
+    bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells);
+    ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells);
     ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells);
     // ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
     double oildistr[2] = { 0.0, 1.0 };
-    ok = ok && append_well_controls(SURFACE_RATE, 8.64297e-05, oildistr, 1, wells);
+    ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells);
     if (!ok) {
         THROW("Something went wrong with well init.");
     }

opm/autodiff/BlackoilPropsAd.hpp

@@ -20,6 +20,7 @@
 #ifndef OPM_BLACKOILPROPSAD_HEADER_INCLUDED
 #define OPM_BLACKOILPROPSAD_HEADER_INCLUDED
 
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 
@@ -37,7 +38,7 @@ namespace Opm
     /// taking an AD type and returning the same. Derivatives are not
     /// returned separately by any method, only implicitly with the AD
     /// version of the methods.
-    class BlackoilPropsAd
+    class BlackoilPropsAd : public BlackoilPropsAdInterface
     {
     public:
         /// Constructor wrapping an opm-core black oil interface.

opm/autodiff/BlackoilPropsAdInterface.cpp

@@ -0,0 +1,24 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+
+Opm::BlackoilPropsAdInterface::~BlackoilPropsAdInterface()
+{
+}

opm/autodiff/BlackoilPropsAdInterface.hpp

@@ -0,0 +1,247 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED
+#define OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+namespace Opm
+{
+
+    /// This class is intended to present a fluid interface for
+    /// three-phase black-oil that is easy to use with the AD-using
+    /// simulators.
+    ///
+    /// Most methods are available in two overloaded versions, one
+    /// taking a constant vector and returning the same, and one
+    /// taking an AD type and returning the same. Derivatives are not
+    /// returned separately by any method, only implicitly with the AD
+    /// version of the methods.
+    class BlackoilPropsAdInterface
+    {
+    public:
+        /// Virtual destructor for inheritance.
+        virtual ~BlackoilPropsAdInterface();
+
+        ////////////////////////////
+        //      Rock interface    //
+        ////////////////////////////
+
+        /// \return   D, the number of spatial dimensions.
+        virtual int numDimensions() const = 0;
+
+        /// \return   N, the number of cells.
+        virtual int numCells() const = 0;
+
+        /// \return   Array of N porosity values.
+        virtual const double* porosity() const = 0;
+
+        /// \return   Array of ND^2 permeability values.
+        ///           The D^2 permeability values for a cell are organized as a matrix,
+        ///           which is symmetric (so ordering does not matter).
+        virtual const double* permeability() const = 0;
+
+
+        ////////////////////////////
+        //      Fluid interface   //
+        ////////////////////////////
+
+        typedef AutoDiff::ForwardBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef std::vector<int> Cells;
+
+
+        // ------ Canonical named indices for each phase ------
+
+        /// Canonical named indices for each phase.
+        enum PhaseIndex { Water = 0, Oil = 1, Gas = 2 };
+
+
+        // ------ Density ------
+
+        /// Densities of stock components at surface conditions.
+        /// \return Array of 3 density values.
+        virtual const double* surfaceDensity() const = 0;
+
+
+        // ------ Viscosity ------
+
+        /// Water viscosity.
+        /// \param[in]  pw     Array of n water pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        V muWat(const V& pw,
+                const Cells& cells) const = 0;
+
+        /// Oil viscosity.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        V muOil(const V& po,
+                const V& rs,
+                const Cells& cells) const = 0;
+
+        /// Gas viscosity.
+        /// \param[in]  pg     Array of n gas pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        V muGas(const V& pg,
+                const Cells& cells) const = 0;
+
+        /// Water viscosity.
+        /// \param[in]  pw     Array of n water pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        ADB muWat(const ADB& pw,
+                  const Cells& cells) const = 0;
+
+        /// Oil viscosity.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        ADB muOil(const ADB& po,
+                  const ADB& rs,
+                  const Cells& cells) const = 0;
+
+        /// Gas viscosity.
+        /// \param[in]  pg     Array of n gas pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n viscosity values.
+        virtual
+        ADB muGas(const ADB& pg,
+                  const Cells& cells) const = 0;
+
+
+        // ------ Formation volume factor (b) ------
+
+        /// Water formation volume factor.
+        /// \param[in]  pw     Array of n water pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        V bWat(const V& pw,
+               const Cells& cells) const = 0;
+
+        /// Oil formation volume factor.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        V bOil(const V& po,
+               const V& rs,
+               const Cells& cells) const = 0;
+
+        /// Gas formation volume factor.
+        /// \param[in]  pg     Array of n gas pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        V bGas(const V& pg,
+               const Cells& cells) const = 0;
+
+        /// Water formation volume factor.
+        /// \param[in]  pw     Array of n water pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        ADB bWat(const ADB& pw,
+                 const Cells& cells) const = 0;
+
+        /// Oil formation volume factor.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        ADB bOil(const ADB& po,
+                 const ADB& rs,
+                 const Cells& cells) const = 0;
+
+        /// Gas formation volume factor.
+        /// \param[in]  pg     Array of n gas pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n formation volume factor values.
+        virtual
+        ADB bGas(const ADB& pg,
+                 const Cells& cells) const = 0;
+
+
+        // ------ Rs bubble point curve ------
+#if 0
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        V rsMax(const V& po,
+                const Cells& cells) const = 0;
+
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        ADB rsMax(const ADB& po,
+                  const Cells& cells) const = 0;
+#endif
+
+        // ------ Relative permeability ------
+
+        /// Relative permeabilities for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  sg     Array of n gas saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 3 elements, each an array of n relperm values,
+        ///                    containing krw, kro, krg. Use PhaseIndex for indexing into the result.
+        virtual
+        std::vector<V> relperm(const V& sw,
+                               const V& so,
+                               const V& sg,
+                               const Cells& cells) const = 0;
+
+        /// Relative permeabilities for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  sg     Array of n gas saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 3 elements, each an array of n relperm values,
+        ///                    containing krw, kro, krg. Use PhaseIndex for indexing into the result.
+        virtual
+        std::vector<ADB> relperm(const ADB& sw,
+                                 const ADB& so,
+                                 const ADB& sg,
+                                 const Cells& cells) const = 0;
+
+    };
+
+} // namespace Opm
+
+#endif // OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED

opm/autodiff/ImpesTPFAAD.hpp

@@ -23,6 +23,7 @@
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
@@ -102,11 +103,11 @@ namespace {
 namespace Opm {
 
 
-    template <class BOFluid, class GeoProps>
+    template <class GeoProps>
     class ImpesTPFAAD {
     public:
         ImpesTPFAAD(const UnstructuredGrid& grid ,
-                    const BOFluid&          fluid,
+                    const BlackoilPropsAdInterface&          fluid,
                     const GeoProps&         geo  ,
                     const Wells&            wells,
                     const LinearSolverInterface& linsolver)
@@ -187,8 +188,6 @@ namespace Opm {
         ImpesTPFAAD(const ImpesTPFAAD& rhs);
         ImpesTPFAAD& operator=(const ImpesTPFAAD& rhs);
 
-        // typedef PressureDependentFluidData<double, BOFluid> PDepFData;
-        // typedef typename PDepFData::ADB ADB;
         typedef AutoDiff::ForwardBlock<double> ADB;
         typedef ADB::V V;
         typedef ADB::M M;
@@ -198,7 +197,7 @@ namespace Opm {
                              Eigen::RowMajor> DataBlock;
 
         const UnstructuredGrid&      grid_;
-        const BOFluid&               fluid_;
+        const BlackoilPropsAdInterface&               fluid_;
         const GeoProps&              geo_ ;
         const Wells&                 wells_;
         const LinearSolverInterface& linsolver_;
@@ -211,9 +210,9 @@ namespace Opm {
         std::vector<V>               kr_;
         std::vector<V>               well_kr_;
 
-        enum { Water = BOFluid::Water,
+        enum { Water = BlackoilPropsAdInterface::Water,
-               Oil = BOFluid::Oil,
+               Oil = BlackoilPropsAdInterface::Oil,
-               Gas = BOFluid::Gas };
+               Gas = BlackoilPropsAdInterface::Gas };
 
 
         void