Work in progress on ensuring units are done correctly.

dune/porsol/blackoil/fluid/BlackoilPVT.cpp

@@ -38,30 +38,25 @@ namespace Opm
         typedef std::vector<std::vector<std::vector<double> > > table_t;
 	region_number_ = 0;
 
-
-
-
-
-
-
 	// Surface densities. Accounting for different orders in eclipse and our code.
 	if (parser.hasField("DENSITY")) {
 	    const int region_number = 0;
 	    enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
-	    const std::vector<std::vector<double> > d_tmp = 
+	    const std::vector<double>& d = parser.getDENSITY().densities_[region_number];
-		parser.getDENSITY().densities_;
+            const double du = parser.units().density;
-	    densities_[Aqua] = d_tmp[region_number][ECL_water];
+            using namespace Dune::unit;
-	    densities_[Liquid] = d_tmp[region_number][ECL_oil];
+	    densities_[Aqua] = convert::from(d[ECL_water], du);
-	    densities_[Vapour] = d_tmp[region_number][ECL_gas];
+	    densities_[Vapour] = convert::from(d[ECL_gas], du);
+	    densities_[Liquid] = convert::from(d[ECL_oil], du);
 	} else {
 	    THROW("Input is missing DENSITY\n");
 	}
 
         // Water PVT
         if (parser.hasField("PVTW")) {
-            water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_));
+            water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_, parser.units()));
         } else {
-            water_props_.reset(new MiscibilityWater(3e-4)); // Default is 0.3 cP.
+            water_props_.reset(new MiscibilityWater(0.5*Dune::prefix::centi*Dune::unit::Poise)); // Eclipse 100 default 
         }
 
         // Oil PVT

dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp

@@ -233,13 +233,11 @@ public:
     { return temperature_; };
 
 public:
-    Scalar fugacity_[numComponents];
-    Scalar moleFrac_[numPhases][numComponents];
-    Scalar phaseConcentration_[numPhases];
-    Scalar meanMolarMass_[numPhases];
-    Scalar phasePressure_[numPhases];
-    Scalar saturation_[numPhases];
     Scalar temperature_;
+    Scalar surface_volume_[numComponents];
+    Scalar phase_pressure_[numPhases];
+    Scalar phase_volume_[numPhases];
+    Scalar saturation_[numPhases];
 };
 
 } // end namespace Opm

dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp

@@ -1,3 +1,17 @@
+/*****************************************************************************
+ *   Copyright (C) 2009 by Andreas Lauser
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version, as long as this copyright notice    *
+ *   is included in its original form.                                       *
+ *                                                                           *
+ *   This program is distributed WITHOUT ANY WARRANTY.                       *
+ *****************************************************************************/
 /*
   Copyright 2010 SINTEF ICT, Applied Mathematics.
 
@@ -28,11 +42,11 @@ namespace Opm
 {
 
 // Forward declaration needed by associated parameters classes.
-template <class ScalarT, class ParamsT>
+template <class ParamT>
-class FluidMatrixInteractionBlackoil;
+class FluidSystemBlackoil;
 
-
+// Parameters for the black oil fluid system.
-class FluidSystemBlackoilParametersNonmiscible
+class FluidSystemBlackoilParameters
 {
 public:
     void init(const Dune::EclipseGridParser& ep)
@@ -40,6 +54,8 @@ public:
         pvt_.init(ep);
     }
 private:
+    template <class P>
+    friend class FluidSystemBlackoil;
     BlackoilPVT pvt_;
 };
 
@@ -47,7 +63,7 @@ private:
 /*!
  * \brief A black oil fluid system.
  */
-template <class ParamsT = FluidSystemBlackoilParametersNonmiscible>
+template <class ParamsT = FluidSystemBlackoilParameters>
 class FluidSystemBlackoil
 {
 public:
@@ -60,6 +76,7 @@ public:
     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.
@@ -96,15 +113,15 @@ public:
     }
 
 
-#if 0
     /*!
      * \brief Return the molar mass of a component [kg/mol].
      */
     static Scalar molarMass(int compIdx)
     {
-        
+        throw std::logic_error("molarMass() not implemented.");
     }
 
+
     /*!
      * \brief Given a phase's composition, temperature, pressure, and
      *        the partial pressures of all components, return its
@@ -116,31 +133,7 @@ public:
                                Scalar pressure,
                                const FluidState &fluidState)
     {
-        if (phaseIdx == lPhaseIdx) {
+        throw std::logic_error("phaseDensity() not implemented.");
-            // See: Ochs 2008
-            // \todo: proper citation
-            Scalar rholH2O = H2O::liquidDensity(temperature, pressure);
-            Scalar clH2O = rholH2O/H2O::molarMass();
-
-            // this assumes each nitrogen molecule displaces exactly one
-            // water molecule in the liquid
-            return
-                clH2O*(H2O::molarMass()*fluidState.moleFrac(lPhaseIdx, H2OIdx)
-                       +
-                       N2::molarMass()*fluidState.moleFrac(lPhaseIdx, N2Idx));
-        }
-        else if (phaseIdx == gPhaseIdx) {
-            Scalar fugH2O = 
-                fluidState.moleFrac(gPhaseIdx, H2OIdx)  *
-                fluidState.phasePressure(gPhaseIdx);
-            Scalar fugN2 = 
-                fluidState.moleFrac(gPhaseIdx, N2Idx)  *
-                fluidState.phasePressure(gPhaseIdx);
-            return
-                H2O::gasDensity(temperature, fugH2O) +
-                N2::gasDensity(temperature, fugN2);
-        }
-        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
     }
 
     /*!
@@ -153,151 +146,41 @@ public:
                                  Scalar pressure,
                                  const FluidState &fluidState)
     {
-        if (phaseIdx == lPhaseIdx) {
+        throw std::logic_error("phaseViscosity() not implemented.");
-            // assume pure water for the liquid phase
-            // TODO: viscosity of mixture
-            return H2O::liquidViscosity(temperature,
-                                        pressure);
-        }
-        else {
-            return N2::gasViscosity(temperature,
-                                    pressure);
-            /* Wilke method. See:
-             *
-             * S.O.Ochs: "Development of a multiphase multicomponent
-             * model for PEMFC - Technical report: IRTG-NUPUS",
-             * University of Stuttgart, 2008
-             *
-             * and:
-             *
-             * See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
-             * edition, McGraw-Hill, 1987, 407-410
-             */
-            Scalar muResult = 0;
-            const Scalar mu[numComponents] = {
-                H2O::gasViscosity(temperature,
-                                  H2O::vaporPressure(temperature)),
-                N2::gasViscosity(temperature,
-                                 pressure)
-            };
-            // molar masses
-            const Scalar M[numComponents] = {
-                H2O::molarMass(),
-                N2::molarMass()
-            };
-
-            for (int i = 0; i < numComponents; ++i) {
-                Scalar divisor = 0;
-                for (int j = 0; j < numComponents; ++j) {
-                    Scalar phiIJ = 1 + sqrt(mu[i]/mu[j] *
-                                            pow(M[i]/M[j], 1/4.0));
-                    phiIJ *= phiIJ;
-                    phiIJ /= sqrt(8*(1 + M[i]/M[j]));
-                    divisor += fluidState.moleFrac(phaseIdx, j)*phiIJ;
-                }
-                muResult += fluidState.moleFrac(phaseIdx, i)*mu[i] / divisor;
-            }
-
-            return muResult;
-        }
     }
 
 
     /*!
-     * \brief Assuming the composition of a single phase and the
+     * \brief Assuming the surface volumes and the pressure of all
-     *        pressure of all phases is known or that all phases are
+     *        phases are known, compute the phase volumes and
-     *        present, compute the thermodynamic equilibrium from the
+     *        saturations.
-     *        temperature and phase pressures. If the known phase
-     *        index
-     *
      */
     template <class FluidState>
-    static void computeEquilibrium(FluidState &fluidState,
+    static void computeEquilibrium(FluidState& fluid_state)
-                                   int knownPhaseIdx = -1)
     {
-        const Scalar T = fluidState.temperature();
+        // Get B and R factors.
-        const Scalar pg = fluidState.phasePressure(gPhaseIdx);
+        const double* p = fluid_state.phase_pressure_;
-        const Scalar pl = fluidState.phasePressure(lPhaseIdx);
+        const double* 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);
+        B[Liquid] = params().pvt_.B(p[Liquid], z, Liquid);
+        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);
 
-        const Scalar betaH2O = H2O::vaporPressure(T);
+        // Update phase volumes.
-        const Scalar betaN2 = BinaryCoeff::H2O_N2::henry(T);
+        double* 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;
+        u[Liquid] = B[Aqua]*(z[Oil] - R[Vapour]*z[Gas])/detR;
 
-        if (knownPhaseIdx < 0)
+        // Update saturations.
-        {
+        double sumu = u[Aqua] + u[Vapour] + u[Liquid];
-            // we only have all phase pressures and temperature and
+        double* s = fluid_state.saturation_;
-            // know that all phases are present
+        for (int i = 0; i < 3; ++i) {
-            Scalar xlH2O = (pg - betaN2)/(betaH2O - betaN2);
+            s[i] = u[i]/sumu;
-            Scalar xlN2 = 1 - xlH2O;
-            Scalar rhol = liquidPhaseDensity_(T, pl, xlH2O, xlN2);
-
-            Scalar cgH2O = H2O::gasDensity(T, betaH2O*xlH2O)/H2O::molarMass();
-            Scalar cgN2 = N2::gasDensity(T, betaN2*xlN2)/N2::molarMass();
-
-            Scalar xgH2O = cgH2O/(cgH2O + cgN2);
-            Scalar xgN2 = cgN2/(cgH2O + cgN2);
-
-            // set the liquid phase composition
-            SettablePhase liquid;
-            liquid.moleFrac_[H2OIdx] = xlH2O;
-            liquid.moleFrac_[N2Idx] = xlN2;
-            liquid.pressure_ = pl;
-            liquid.density_ = rhol;
-            liquid.xToX(); // compute mass fractions from mole fractions
-            fluidState.assignPhase(lPhaseIdx, liquid);
-
-            // set the gas phase composition
-            SettablePhase gas;
-            gas.moleFrac_[H2OIdx] = xgH2O;
-            gas.moleFrac_[N2Idx] = xgN2;
-            gas.pressure_ = pg;
-            gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass();
-            gas.xToX(); // compute mass fractions from mole fractions
-            fluidState.assignPhase(gPhaseIdx, gas);
-        }
-        else if (knownPhaseIdx == lPhaseIdx) {
-            // the composition of the liquid phase is given
-
-            // retrieve the known mole fractions from the fluid state
-            Scalar xlH2O = fluidState.moleFrac(lPhaseIdx, H2OIdx);
-            Scalar xlN2 = fluidState.moleFrac(lPhaseIdx, N2Idx);
-
-            // calculate the component contentrations in the gas phase
-            Scalar pH2O = betaH2O*xlH2O; // fugacity of water
-            Scalar pN2 = betaN2*xlN2; // fugacity of nitrogen
-            Scalar cgH2O = H2O::gasDensity(T, pH2O)/H2O::molarMass();
-            Scalar cgN2 = N2::gasDensity(T, pN2)/N2::molarMass();
-
-            // convert concentrations to mole fractions
-            Scalar xgH2O = cgH2O/(cgH2O + cgN2) * (pH2O + pN2)/pg;
-            Scalar xgN2 = cgN2/(cgH2O + cgN2) * (pH2O + pN2)/pg;
-
-            // set gas phase composition
-            SettablePhase gas;
-            gas.moleFrac_[H2OIdx] = xgH2O;
-            gas.moleFrac_[N2Idx] = xgN2;
-            gas.pressure_ = pg;
-            gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass();
-            gas.xToX(); // update mass fractions from mole fractions
-            fluidState.assignPhase(gPhaseIdx, gas);
-        }
-        else if (knownPhaseIdx == gPhaseIdx) {
-            // the composition of the gas phase is given
-
-            Scalar xgH2O = fluidState.moleFrac(gPhaseIdx, H2OIdx);
-            Scalar xgN2 = fluidState.moleFrac(gPhaseIdx, N2Idx);
-            Scalar pgH2O = pg*xgH2O;
-            Scalar pgN2 = pg*xgN2;
-
-            Scalar xlH2O = pgH2O/betaH2O;
-            Scalar xlN2 = pgN2/betaN2;
-
-            SettablePhase liquid;
-            liquid.moleFrac_[H2OIdx] = xlH2O;
-            liquid.moleFrac_[N2Idx] = xlN2;
-            liquid.pressure_ = pl;
-            liquid.density_ = liquidPhaseDensity_(T, pl, xlH2O, xlN2);
-            liquid.xToX(); // update mass fractions from mole fractions
-            fluidState.assignPhase(lPhaseIdx, liquid);
         }
     }
 
@@ -326,28 +209,8 @@ public:
                                 Scalar pressure,
                                 const FluidState &state)
     {
-        if (phaseIdx == gPhaseIdx) {
+        throw std::logic_error("activityCoeff() not implemented.");
-            return pressure;
+        return 0.0;
-            Scalar fugH2O = std::max(1e-3, state.fugacity(H2OIdx));
-            Scalar fugN2 = std::max(1e-3, state.fugacity(N2Idx));
-            Scalar cH2O = H2O::gasDensity(temperature, fugH2O) / H2O::molarMass();
-            Scalar cN2 = N2::gasDensity(temperature, fugN2) / N2::molarMass();
-
-            Scalar alpha = (fugH2O + fugN2)/pressure;
-
-            if (compIdx == H2OIdx)
-                return fugH2O/(alpha*cH2O/(cH2O + cN2));
-            else if (compIdx == N2Idx)
-                return fugN2/(alpha*cN2/(cH2O + cN2));
-
-            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
-        }
-
-        switch (compIdx) {
-        case H2OIdx: return H2O::vaporPressure(temperature);
-        case N2Idx: return BinaryCoeff::H2O_N2::henry(temperature);
-        };
-        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
     }
 
     /*!
@@ -363,49 +226,8 @@ public:
                             Scalar pressure,
                             const FluidState &fluidState)
     {
-        if (compIIdx > compJIdx)
+        throw std::logic_error("diffCoeff() not implemented.");
-            std::swap(compIIdx, compJIdx);
+        return 0.0;
-
-#ifndef NDEBUG
-        if (compIIdx == compJIdx ||
-            phaseIdx > numPhases - 1 ||
-            compJIdx > numComponents - 1)
-        {
-            DUNE_THROW(Dune::InvalidStateException,
-                       "Binary diffusion coefficient of components "
-                       << compIIdx << " and " << compJIdx
-                       << " in phase " << phaseIdx << " is undefined!\n");
-        }
-#endif
-
-
-        switch (phaseIdx) {
-        case lPhaseIdx:
-            switch (compIIdx) {
-            case H2OIdx:
-                switch (compJIdx) {
-                case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(temperature,
-                                                                        pressure);
-                }
-            default:
-                DUNE_THROW(Dune::InvalidStateException,
-                           "Binary diffusion coefficients of trace "
-                           "substances in liquid phase is undefined!\n");
-            }
-        case gPhaseIdx:
-            switch (compIIdx) {
-            case H2OIdx:
-                switch (compJIdx) {
-                case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(temperature,
-                                                                     pressure);
-                }
-            }
-        }
-
-        DUNE_THROW(Dune::InvalidStateException,
-                   "Binary diffusion coefficient of components "
-                   << compIIdx << " and " << compJIdx
-                   << " in phase " << phaseIdx << " is undefined!\n");
     };
 
     /*!
@@ -418,35 +240,8 @@ public:
                                 Scalar pressure,
                                 const FluidState &fluidState)
     {
-        if (phaseIdx == lPhaseIdx) {
+        throw std::logic_error("phaseEnthalpy() not implemented.");
-            Scalar cN2 = fluidState.concentration(lPhaseIdx, N2Idx);
+        return 0.0;
-            Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass());
-
-            // TODO: correct way to deal with the solutes???
-            return
-                fluidState.massFrac(lPhaseIdx, H2OIdx)*
-                H2O::liquidEnthalpy(temperature, pressure)
-                +
-                fluidState.massFrac(lPhaseIdx, N2Idx)*
-                N2::gasEnthalpy(temperature, pN2);
-        }
-        else {
-            Scalar cH2O = fluidState.concentration(gPhaseIdx, H2OIdx);
-            Scalar cN2 = fluidState.concentration(gPhaseIdx, N2Idx);
-            
-            Scalar pH2O = H2O::gasPressure(temperature, cH2O*H2O::molarMass());
-            Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass());
-
-            Scalar result = 0;
-            result +=
-                H2O::gasEnthalpy(temperature, pH2O) *
-                fluidState.massFrac(gPhaseIdx, H2OIdx);
-            result +=
-                N2::gasEnthalpy(temperature, pN2) *
-                fluidState.massFrac(gPhaseIdx, N2Idx);
-
-            return result;
-        }
     }
 
     /*!
@@ -459,33 +254,11 @@ public:
                                       Scalar pressure,
                                       const FluidState &fluidState)
     {
-        return
+        throw std::logic_error("phaseInternalEnergy() not implemented.");
-            phaseEnthalpy(phaseIdx, temperature, pressure, fluidState) -
+        return 0.0;
-            pressure/phaseDensity(phaseIdx, temperature, pressure, fluidState);
     }
 
 private:
-    static Scalar liquidPhaseDensity_(Scalar T, Scalar pl, Scalar xlH2O, Scalar xlN2)
-    {
-        // See: Ochs 2008
-        // \todo: proper citation
-        Scalar rholH2O = H2O::liquidDensity(T, pl);
-        Scalar clH2O = rholH2O/H2O::molarMass();
-
-        // this assumes each nitrogen molecule displaces exactly one
-        // water molecule in the liquid
-        return
-            clH2O*(xlH2O*H2O::molarMass()
-                   +
-                   xlN2*N2::molarMass());
-    }
-
-    static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
-    {
-        return H2O::gasDensity(T, pg*xgH2O) + N2::gasDensity(T, pg*xgN2);
-    };
-
-#endif
 
     static Params& params()
     {

dune/porsol/blackoil/fluid/MiscibilityWater.hpp

@@ -36,6 +36,7 @@
 
 #include "MiscibilityProps.hpp"
 #include <dune/common/ErrorMacros.hpp>
+#include <dune/common/Units.hpp>
 
 // Forward declaration.
 class PVTW;
@@ -46,24 +47,26 @@ namespace Opm
     {
     public:
         typedef std::vector<std::vector<double> > table_t;
-	MiscibilityWater(const table_t& pvtw)
+	MiscibilityWater(const table_t& pvtw, const Dune::EclipseUnits& units)
         {
 	    const int region_number = 0;
 	    if (pvtw.size() != 1) {
 		THROW("More than one PVD-region");
 	    }
-            double b = pvtw[region_number][1];
+            using namespace Dune::unit;
-            double comp = pvtw[region_number][2];
+            ref_press_ = convert::from(pvtw[region_number][0], units.pressure);
-            double visc = pvtw[region_number][3];
+            ref_B_ = convert::from(pvtw[region_number][1], units.liqvol_r/units.liqvol_s);
-            const double VISCOSITY_UNIT = 1e-3;
+            comp_ = convert::from(pvtw[region_number][2], units.compressibility);
-            if (b == 1.0 && comp == 0) {
+            viscosity_ = convert::from(pvtw[region_number][3], units.viscosity);
-                viscosity_ = visc*VISCOSITY_UNIT;
+            if (pvtw[region_number].size() > 4 && pvtw[region_number][4] != 0.0) {
-            } else {
+                THROW("MiscibilityWater does not support 'viscosibility'.");
-                THROW("Compressible water not implemented.");
             }
         }
 	MiscibilityWater(double visc)
-            : viscosity_(visc)
+            : ref_press_(0.0),
+              ref_B_(1.0),
+              comp_(0.0),
+              viscosity_(visc)
         {
         }
 	virtual ~MiscibilityWater()
@@ -91,6 +94,9 @@ namespace Opm
             return 0.0;
         }
     private:
+        double ref_press_;
+        double ref_B_;
+        double comp_;
         double viscosity_;
     };