A large number of additions to start testing compressible tpfa-solver.

dune/porsol/blackoil/BlackoilFluid.hpp

@@ -0,0 +1,68 @@
+/*
+  Copyright 2010 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_BLACKOILFLUID_HEADER_INCLUDED
+#define OPM_BLACKOILFLUID_HEADER_INCLUDED
+
+
+#include <dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp>
+#include <dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp>
+#include <dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp>
+#include <dune/common/EclipseGridParser.hpp>
+#include <dune/common/fvector.hh>
+
+
+namespace Opm
+{
+
+    class BlackoilFluid
+    {
+    public:
+        enum { numPhases = 3 };
+        enum { numComponents = 3 };
+        typedef Dune::FieldVector<double, numPhases> PhaseVec;
+        typedef Dune::FieldVector<double, numComponents> CompVec;
+
+        void init(const Dune::EclipseGridParser& parser)
+        {
+            fmi_params_.init(parser);
+            FluidSystemBlackoil<>::init(parser);
+        }
+
+        FluidStateBlackoil computeState(PhaseVec phase_pressure, CompVec z)
+        {
+            FluidStateBlackoil state;
+            state.temperature_ = 300;
+            state.phase_pressure_ = phase_pressure;
+            state.surface_volume_ = z;
+            FluidSystemBlackoil<>::computeEquilibrium(state); // Sets everything but relperm and mobility.
+            FluidMatrixInteractionBlackoil<double>::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_);
+            for (int phase = 0; phase < numPhases; ++phase) {
+                state.mobility_[phase] = state.relperm_[phase]/state.viscosity_[phase];
+            }
+            return state;
+        }
+
+    private:
+        FluidMatrixInteractionBlackoilParams<double> fmi_params_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_BLACKOILFLUID_HEADER_INCLUDED

dune/porsol/blackoil/fluid/BlackoilDefs.hpp

@@ -20,6 +20,10 @@
 #ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED
 #define OPM_BLACKOILDEFS_HEADER_INCLUDED
 
+
+#include <dune/common/fvector.hh>
+
+
 namespace Opm
 {
 
@@ -31,6 +35,10 @@ namespace Opm
 
         enum ComponentIndex { Water = 0, Gas = 1, Oil = 2 };
         enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 };
+
+        typedef double Scalar;
+        typedef Dune::FieldVector<Scalar, numComponents> CompVec;
+        typedef Dune::FieldVector<Scalar, numPhases> PhaseVec;
     };
 
 } // namespace Opm

dune/porsol/blackoil/fluid/BlackoilPVT.cpp

@@ -78,32 +78,32 @@ namespace Opm
         }
     }
 
-    BlackoilPVT::surfvol_t BlackoilPVT::surfaceDensities() const
+    BlackoilPVT::CompVec BlackoilPVT::surfaceDensities() const
     {
         return densities_;
     }
 
-    double BlackoilPVT::getViscosity(double press, const surfvol_t& surfvol, PhaseIndex phase) const
+    double BlackoilPVT::getViscosity(double press, const CompVec& surfvol, PhaseIndex phase) const
     {
         return propsForPhase(phase).getViscosity(region_number_, press, surfvol);
     }
 
-    double BlackoilPVT::B(double press, const surfvol_t& surfvol, PhaseIndex phase) const
+    double BlackoilPVT::B(double press, const CompVec& surfvol, PhaseIndex phase) const
     {
         return propsForPhase(phase).B(region_number_, press, surfvol);
     }
 
-    double BlackoilPVT::dBdp(double press, const surfvol_t& surfvol, PhaseIndex phase) const
+    double BlackoilPVT::dBdp(double press, const CompVec& surfvol, PhaseIndex phase) const
     {
         return propsForPhase(phase).dBdp(region_number_, press, surfvol);
     }
 
-    double BlackoilPVT::R(double press, const surfvol_t& surfvol, PhaseIndex phase) const
+    double BlackoilPVT::R(double press, const CompVec& surfvol, PhaseIndex phase) const
     {
         return propsForPhase(phase).R(region_number_, press, surfvol);
     }
 
-    double BlackoilPVT::dRdp(double press, const surfvol_t& surfvol, PhaseIndex phase) const
+    double BlackoilPVT::dRdp(double press, const CompVec& surfvol, PhaseIndex phase) const
     {
         return propsForPhase(phase).dRdp(region_number_, press, surfvol);
     }

dune/porsol/blackoil/fluid/BlackoilPVT.hpp

@@ -33,20 +33,23 @@ namespace Opm
     class BlackoilPVT : public BlackoilDefs
     {
     public:
-        typedef MiscibilityProps::surfvol_t surfvol_t;
-
 	void init(const Dune::EclipseGridParser& ep);
 
-        double getViscosity(double press, const surfvol_t& surfvol,
+        double getViscosity(double press,
+                            const CompVec& surfvol,
 			    PhaseIndex phase) const;
-	surfvol_t surfaceDensities() const;
+	CompVec surfaceDensities() const;
-        double B   (double press, const surfvol_t& surfvol,
+        double B   (double press,
+                    const CompVec& surfvol,
 		    PhaseIndex phase) const;
-        double dBdp(double press, const surfvol_t& surfvol,
+        double dBdp(double press,
+                    const CompVec& surfvol,
 		    PhaseIndex phase) const;
-        double R   (double press, const surfvol_t& surfvol,
+        double R   (double press,
+                    const CompVec& surfvol,
 		    PhaseIndex phase) const;
-        double dRdp(double press, const surfvol_t& surfvol,
+        double dRdp(double press,
+                    const CompVec& surfvol,
 		    PhaseIndex phase) const;
 
 
@@ -57,7 +60,7 @@ namespace Opm
 	boost::scoped_ptr<MiscibilityProps> water_props_;
 	boost::scoped_ptr<MiscibilityProps> oil_props_;
 	boost::scoped_ptr<MiscibilityProps> gas_props_;
-	surfvol_t densities_;
+	CompVec densities_;
     };
 
 }

dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp

@@ -21,53 +21,25 @@
 #define OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
 
 
-#include "FluidSystemBlackoil.hpp"
+#include "BlackoilDefs.hpp"
 
 
 namespace Opm
 {
     /*!
-     * \brief Calculates the phase state from the primary variables in the
+     * \brief Fluid states for a black oil model.
-     *        blackoil model.
      */
     struct FluidStateBlackoil : public BlackoilDefs
     {
-        typedef FluidSystemBlackoil FluidSystem;
-        typedef double Scalar;
-        typedef std::tr1::array<Scalar, numComponents + numPhases> PrimaryVariables; // Surface volumes and phase pressures.
-        typedef FluidMatrixInteractionBlackoil<Scalar> MaterialLaw;
-        typedef typename MaterialLaw::Params MaterialLawParams;
-
         Scalar temperature_;
-        Scalar surface_volume_[numComponents];
+        CompVec surface_volume_;
-        Scalar phase_pressure_[numPhases];
+        PhaseVec phase_pressure_;
-        Scalar phase_volume_[numPhases];
+        PhaseVec phase_volume_;
-        Scalar saturation_[numPhases];
+        Scalar phase_to_comp_[numPhases*numComponents];
-
+        PhaseVec saturation_;
-
+        PhaseVec viscosity_;
-        /*!
+        PhaseVec relperm_;
-         * \brief Update the phase state from the primary variables.
+        PhaseVec mobility_;
-         */
-        void update(const PrimaryVariables& primary_vars,
-                    const MaterialLawParams& material_params,
-                    Scalar temperature)
-        {
-            // Set the temperature.
-            temperature_ = temperature;
-
-            // Set the surface volumes.
-            for (int i = 0; i < numComponents; ++i) {
-                surface_volume_[i] = primary_vars[i];
-            }
-
-            // Set the phase pressures.
-            for (int i = 0; i < numPhases; ++i) {
-                phase_pressure_[i] = primary_vars[numComponents + i];
-            }
-
-            // Compute phase volumes by the fluid system rules.
-            FluidSystem::computeEquilibrium(*this);
-        }
     };
 
 } // end namespace Opm

dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp

@@ -35,6 +35,8 @@
 #define OPM_FLUIDSYSTEMBLACKOIL_HEADER_INCLUDED
 
 #include "BlackoilPVT.hpp"
+#include "BlackoilDefs.hpp"
+#include <dune/porsol/common/Matrix.hpp>
 #include <dune/common/EclipseGridParser.hpp>
 #include <stdexcept>
 
@@ -64,19 +66,10 @@ private:
  * \brief A black oil fluid system.
  */
 template <class ParamsT = FluidSystemBlackoilParameters>
-class FluidSystemBlackoil
+class FluidSystemBlackoil : public BlackoilDefs
 {
 public:
     typedef ParamsT Params;
-    typedef double Scalar;
-
-    enum { numComponents = 3 };
-    enum { numPhases = 3 };
-
-    enum ComponentIndex { Water = 0, Gas = 1, Oil = 2 };
-    enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 };
-
-    typedef BlackoilPVT::surfvol_t FluidVec;
 
     /*!
      * \brief Initialize system from input.
@@ -151,16 +144,16 @@ public:
 
 
     /*!
-     * \brief Assuming the surface volumes and the pressure of all
+     * \brief Assuming the surface volumes and the pressures of all
-     *        phases are known, compute the phase volumes and
+     *        phases are known, compute everything except relperm and
-     *        saturations.
+     *        mobility.
      */
     template <class FluidState>
     static void computeEquilibrium(FluidState& fluid_state)
     {
         // Get B and R factors.
-        const double* p = fluid_state.phase_pressure_;
+        const PhaseVec& p = fluid_state.phase_pressure_;
-        const double* z = fluid_state.surface_volume_;
+        const CompVec& z = fluid_state.surface_volume_;
         double B[3];
         B[Aqua]   = params().pvt_.B(p[Aqua],   z, Aqua);
         B[Vapour] = params().pvt_.B(p[Vapour], z, Vapour);
@@ -168,9 +161,17 @@ public:
         double R[3]; // Only using 2 of them, though.
         R[Vapour] = params().pvt_.B(p[Vapour], z, Vapour);
         R[Liquid] = params().pvt_.B(p[Liquid], z, Liquid);
+        // Set the A matrix (A = RB^{-1})
+        Dune::SharedFortranMatrix A(numComponents, numPhases, fluid_state.phase_to_comp_);
+        zero(A);
+        A(Water, Aqua) = 1.0/B[Aqua];
+        A(Gas, Vapour) = 1.0/B[Vapour];
+        A(Gas, Liquid) = R[Liquid]/B[Liquid];
+        A(Oil, Vapour) = R[Vapour]/B[Vapour];
+        A(Oil, Liquid) = 1.0/B[Liquid];
 
-        // Update phase volumes.
+        // Update phase volumes. This is the same as multiplying with A^{-1}
-        double* u = fluid_state.phase_volume_;
+        PhaseVec& u = fluid_state.phase_volume_;
         double detR = 1.0 - R[Vapour]*R[Liquid];
         u[Aqua] = B[Aqua]*z[Water];
         u[Vapour] = B[Vapour]*(z[Gas] - R[Liquid]*z[Oil])/detR;
@@ -178,10 +179,16 @@ public:
 
         // Update saturations.
         double sumu = u[Aqua] + u[Vapour] + u[Liquid];
-        double* s = fluid_state.saturation_;
+        PhaseVec& s = fluid_state.saturation_;
         for (int i = 0; i < 3; ++i) {
             s[i] = u[i]/sumu;
         }
+
+        // Compute viscosities.
+        PhaseVec& mu = fluid_state.viscosity_;
+        mu[Aqua]   = params().pvt_.getViscosity(p[Aqua],   z, Aqua);
+        mu[Vapour] = params().pvt_.getViscosity(p[Vapour], z, Vapour);
+        mu[Liquid] = params().pvt_.getViscosity(p[Liquid], z, Liquid);
     }
 
     /*!

dune/porsol/blackoil/fluid/MiscibilityProps.hpp

@@ -46,7 +46,7 @@ namespace Opm
     class MiscibilityProps : public BlackoilDefs
     {
     public:
-	typedef std::tr1::array<double, 3> surfvol_t;
+	typedef CompVec surfvol_t;
 
 	MiscibilityProps();
 	virtual ~MiscibilityProps();