fix immiscible flash calculation, add test for it

also add the peng-robinson test to the CMake build

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7595 2fb0f335-1f38-0410-981e-8018bf24f1b0

configure.ac

@@ -94,10 +94,11 @@ AC_CONFIG_FILES([dumux.pc
     test/freeflow/stokes2c/Makefile
     test/freeflow/stokes2cni/Makefile
     test/material/Makefile
-    test/material/tabulation/Makefile
+    test/material/fluidsystems/Makefile
+    test/material/immiscibleflash/Makefile
     test/material/ncpflash/Makefile
     test/material/pengrobinson/Makefile
-    test/material/fluidsystems/Makefile
+    test/material/tabulation/Makefile
     tutorial/Makefile
 ])
 

dumux/material/constraintsolvers/immiscibleflash.hh

@@ -49,7 +49,7 @@ namespace Dumux {
  * This sums up to 2*N unknowns. On the equations side of things, we
  * have:
  *
- * - N total masses
+ * - N total component molarities
  * - 1 The sum of all saturations is 1
  * - N-1 Relations from capillary pressure
  *
@@ -94,29 +94,14 @@ public:
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
             sumMoles += globalMolarities[compIdx];
 
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-            // composition
-            for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
-                fluidState.setMoleFraction(phaseIdx,
-                                           compIdx,
-                                           globalMolarities[compIdx]/sumMoles);
-            
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {          
             // pressure. use atmospheric pressure as initial guess
             fluidState.setPressure(phaseIdx, 2e5);
 
             // saturation. assume all fluids to be equally distributed
             fluidState.setSaturation(phaseIdx, 1.0/numPhases);
         }
-        
-        // set the fugacity coefficients of all components in all phases
-        paramCache.updateAll(fluidState);
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-            for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
-                Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
-                fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
-            }
-        }
-    };
+    }
 
     /*!
      * \brief Calculates the chemical equilibrium from the component
@@ -211,7 +196,7 @@ public:
                 }
                 std::cout << "\n";
             };
-
+            std::cout << "residual: " << b << "\n";
             std::cout << "deltaX: " << deltaX << "\n";
             std::cout << "---------------\n";
             */
@@ -235,7 +220,7 @@ public:
                    "Flash calculation failed."
                    " {c_alpha^kappa} = {" << globalMolarities << "}, T = " 
                    << fluidState.temperature(/*phaseIdx=*/0));
-    };
+    }
 
 
 protected:
@@ -298,7 +283,7 @@ protected:
 
             // deviate the mole fraction of the i-th component
             Scalar x_i = getQuantity_(fluidState, pvIdx);
-            const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
+            const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
             setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
             assert(getQuantity_(fluidState, pvIdx) == x_i + eps);
             
@@ -306,7 +291,6 @@ protected:
             calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
             tmp -= b;
             tmp /= eps;
-
             // store derivative in jacobian matrix
             for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
                 J[eqIdx][pvIdx] = tmp[eqIdx];
@@ -378,8 +362,11 @@ protected:
         Scalar sumSat = 0.0;
         for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
             sumSat += fluidState.saturation(phaseIdx);
+
+        // make sure that the last saturation does not get out of range [0, 1]
+        sumSat = std::min(1.0, std::max(0.0, sumSat));
         fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
-        
+    
         // update the pressures using the material law (saturations
         // and first pressure are already set because it is implicitly
         // solved for.)
@@ -391,7 +378,7 @@ protected:
                                    + (pC[phaseIdx] - pC[0]));
 
         // update the parameter cache
-        paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature);
+        paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature|ParameterCache::Composition);
 
         // update all densities
         for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
@@ -417,12 +404,13 @@ protected:
             int phaseIdx = 0;
             return fs.pressure(phaseIdx);
         }
-        // first M - 1 saturations
-        else { // if (pvIdx < numPhases)
+        // saturations
+        else { 
+            assert(pvIdx < numPhases);
             int phaseIdx = pvIdx - 1;
             return fs.saturation(phaseIdx);
         }
-    };
+    }
 
     // set a quantity in the fluid state
     template <class MaterialLaw, class FluidState>
@@ -434,7 +422,8 @@ protected:
     {
         assert(0 <= pvIdx && pvIdx < numEq);
 
-        if (pvIdx < 1) { // <- first pressure
+        if (pvIdx < 1) {
+            // -> first pressure
             Scalar delta = value - fs.pressure(0);
 
             // set all pressures. here we assume that the capillary
@@ -449,12 +438,15 @@ protected:
                 fs.setDensity(phaseIdx, rho);
             }
         }
-        else if (pvIdx < numPhases) { // <- first M - 1 saturations
+        else {
+            assert(pvIdx < numPhases);
+
+            // -> first M-1 saturations
             Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
             fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
 
             // set last saturation (-> minus the change of the
-            // satuation of the other phase)
+            // saturation of the other phases)
             fs.setSaturation(/*phaseIdx=*/numPhases - 1, 
                              fs.saturation(numPhases - 1) - delta);
 
@@ -474,10 +466,7 @@ protected:
                 fs.setDensity(phaseIdx, rho);
             }
         }
-        else {
-            assert(false);
-        }
-    };
+    }
 
     // set a quantity in the fluid state
     template <class FluidState>
@@ -490,12 +479,17 @@ protected:
             int phaseIdx = 0;
             fs.setPressure(phaseIdx, value);
         }
-        // first M - 1 saturations
-        else if (pvIdx < numPhases) {
+        // saturations
+        else {
+            assert(pvIdx < numPhases);
             int phaseIdx = pvIdx - 1;
+
+            // make sure that the first M-1 saturations does not get
+            // out of the [0, 1] intervall
+            value = std::min(1.0, std::max(0.0, value));
             fs.setSaturation(phaseIdx, value);
         }
-    };
+    }
 
     template <class FluidState>
     static Scalar quantityWeight_(const FluidState &fs, int pvIdx)
@@ -504,9 +498,11 @@ protected:
         if (pvIdx < 1)
             return 1.0/1e5;
         // first M - 1 saturations
-        else // if (pvIdx < numPhases)
+        else {
+            assert(pvIdx < numPhases);
             return 1.0;
-    };
+        }
+    }
 };
 
 } // end namespace Dumux

dumux/material/fluidsystems/spe5parametercache.hh

@@ -52,7 +52,6 @@ class Spe5ParameterCache
     typedef Dumux::PengRobinson<Scalar> PengRobinson;
 
     enum { numPhases = FluidSystem::numPhases };
-    enum { numComponents = FluidSystem::numComponents };
 
     enum { wPhaseIdx = FluidSystem::wPhaseIdx };
     enum { oPhaseIdx = FluidSystem::oPhaseIdx };

test/material/CMakeLists.txt

@@ -2,6 +2,8 @@
 # (they also don't compile if using the old build system)!
 
 add_subdirectory("fluidsystems")
+add_subdirectory("immiscibleflash")
 add_subdirectory("ncpflash")
+add_subdirectory("pengrobinson")
 add_subdirectory("tabulation")
 

test/material/Makefile.am

-SUBDIRS = fluidsystems ncpflash pengrobinson tabulation 
+SUBDIRS = fluidsystems immiscibleflash ncpflash pengrobinson tabulation 
 
 include $(top_srcdir)/am/global-rules
 

test/material/immiscibleflash/CMakeLists.txt

+# build target for the simple twophase lens problem
+ADD_EXECUTABLE("test_immiscibleflash" test_ncpflash.cc)
+TARGET_LINK_LIBRARIES("test_immiscibleflash" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags 
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+

test/material/immiscibleflash/Makefile.am

+check_PROGRAMS = test_immiscibleflash
+
+test_immiscibleflash_SOURCES = test_immiscibleflash.cc
+
+include $(top_srcdir)/am/global-rules

test/material/immiscibleflash/test_immiscibleflash.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   Copyright (C) 2012 by Andreas Lauser                                    *
+ *   Institute for Modelling Hydraulic and Environmental Systems             *
+ *   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.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This is a program to test the flash calculation which uses
+ *        non-linear complementarity problems (NCP)
+ *
+ * A flash calculation determines the pressures, saturations and
+ * composition of all phases given the total mass (or, as in this case
+ * the total number of moles) in a given amount of pore space.
+ */
+#include "config.h"
+
+#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+#include <dumux/material/constraintsolvers/immiscibleflash.hh>
+
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+
+#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
+
+#include <dumux/material/fluidmatrixinteractions/Mp/Mplinearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/Mp/2padapter.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+template <class Scalar, class FluidState>
+void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
+{
+    enum { numPhases = FluidState::numPhases };
+    enum { numComponents = FluidState::numComponents };
+
+    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+        Scalar error;
+
+        // check the pressures
+        error = 1 - fsRef.pressure(phaseIdx)/fsFlash.pressure(phaseIdx);
+        if (std::abs(error) > 1e-6) {
+            std::cout << "pressure error phase " << phaseIdx << ": "
+                      << fsFlash.pressure(phaseIdx)  << " flash vs "
+                      << fsRef.pressure(phaseIdx) << " reference"
+                      << " error=" << error << "\n";
+        }
+
+        // check the saturations
+        error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx);
+        if (std::abs(error) > 1e-6)
+            std::cout << "saturation error phase " << phaseIdx << ": " 
+                      << fsFlash.saturation(phaseIdx) << " flash vs "
+                      << fsRef.saturation(phaseIdx) << " reference"
+                      << " error=" << error << "\n";
+
+        // check the compositions
+        for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
+            error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx);
+            if (std::abs(error) > 1e-6)
+                std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": " 
+                          << fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs "
+                          << fsRef.moleFraction(phaseIdx, compIdx) << " reference"
+                          << " error=" << error << "\n";
+        }
+    }
+}
+
+template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
+void checkImmiscibleFlash(const FluidState &fsRef, 
+                          typename MaterialLaw::Params &matParams)
+{
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+
+    // calculate the total amount of stuff in the reference fluid
+    // phase
+    ComponentVector globalMolarities(0.0);
+    for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            globalMolarities[compIdx] +=
+                fsRef.saturation(phaseIdx)*fsRef.molarity(phaseIdx, compIdx);
+        }
+    }
+
+    // initialize the fluid state for the flash calculation
+    typedef Dumux::ImmiscibleFlash<Scalar, FluidSystem> ImmiscibleFlash;
+    FluidState fsFlash;
+
+    fsFlash.setTemperature(fsRef.temperature(/*phaseIdx=*/0));
+
+    // run the flash calculation
+    typename FluidSystem::ParameterCache paramCache;
+    ImmiscibleFlash::guessInitial(fsFlash, paramCache, globalMolarities);
+    ImmiscibleFlash::template solve<MaterialLaw>(fsFlash, paramCache, matParams, globalMolarities);
+
+    // compare the "flashed" fluid state with the reference one
+    checkSame<Scalar>(fsRef, fsFlash);
+}
+
+
+template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
+void completeReferenceFluidState(FluidState &fs, 
+                                 typename MaterialLaw::Params &matParams,
+                                 int refPhaseIdx)
+{
+    enum { numPhases = FluidSystem::numPhases };
+
+    typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
+    typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
+    
+    int otherPhaseIdx = 1 - refPhaseIdx;
+
+    // calculate the other saturation
+    fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
+
+    // calulate the capillary pressure
+    PhaseVector pC;
+    MaterialLaw::capillaryPressures(pC, matParams, fs);
+    fs.setPressure(otherPhaseIdx, 
+                   fs.pressure(refPhaseIdx)
+                   + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
+
+    // set all phase densities
+    typename FluidSystem::ParameterCache paramCache;
+    paramCache.updateAll(fs);
+    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+        Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
+        fs.setDensity(phaseIdx, rho);
+    }
+}
+
+
+int main()
+{
+    typedef double Scalar;
+    typedef Dumux::FluidSystems::H2ON2<Scalar> FluidSystem;
+    typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> ImmiscibleFluidState;
+
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+    enum { lPhaseIdx = FluidSystem::lPhaseIdx };
+    enum { gPhaseIdx = FluidSystem::gPhaseIdx };
+
+    enum { H2OIdx = FluidSystem::H2OIdx };
+    enum { N2Idx = FluidSystem::N2Idx };
+
+    typedef Dumux::RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+    typedef Dumux::EffToAbsLaw<EffMaterialLaw> TwoPMaterialLaw;
+    typedef Dumux::TwoPAdapter<lPhaseIdx, TwoPMaterialLaw> MaterialLaw;
+    typedef MaterialLaw::Params MaterialLawParams;
+
+    Scalar T = 273.15 + 25;
+    
+    // initialize the tables of the fluid system
+    Scalar Tmin = T - 1.0;
+    Scalar Tmax = T + 1.0;
+    int nT = 3;
+    
+    Scalar pmin = 0.0;
+    Scalar pmax = 1.25 * 2e6;
+    int np = 100;
+
+    FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
+
+    // set the parameters for the capillary pressure law
+    MaterialLawParams matParams;
+    matParams.setSwr(0.0);
+    matParams.setSnr(0.0);
+    matParams.setPe(0);
+    matParams.setLambda(2.0);
+
+    ImmiscibleFluidState fsRef;
+    
+    // create an fluid state which is consistent
+
+    // set the fluid temperatures
+    fsRef.setTemperature(T);
+    
+    ////////////////
+    // only liquid
+    ////////////////
+    std::cout << "testing single-phase liquid\n";
+
+    // set liquid saturation and pressure
+    fsRef.setSaturation(lPhaseIdx, 1.0);
+    fsRef.setPressure(lPhaseIdx, 1e6);
+    
+    // set the remaining parameters of the reference fluid state
+    completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
+    
+    // check the flash calculation
+    checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+    ////////////////
+    // only gas
+    ////////////////
+    std::cout << "testing single-phase gas\n";
+
+    // set gas saturation and pressure
+    fsRef.setSaturation(gPhaseIdx, 1.0);
+    fsRef.setPressure(gPhaseIdx, 1e6);
+
+    // set the remaining parameters of the reference fluid state
+    completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, gPhaseIdx);
+    
+    // check the flash calculation
+    checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+    ////////////////
+    // both phases
+    ////////////////
+    std::cout << "testing two-phase\n";
+
+    // set liquid saturation and pressure
+    fsRef.setSaturation(lPhaseIdx, 0.5);
+    fsRef.setPressure(lPhaseIdx, 1e6);
+
+    // set the remaining parameters of the reference fluid state
+    completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
+       
+    // check the flash calculation
+    checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+    ////////////////
+    // with capillary pressure
+    ////////////////
+    std::cout << "testing two-phase with capillary pressure\n";
+
+    MaterialLawParams matParams2;
+    matParams2.setSwr(0.0);
+    matParams2.setSnr(0.0);
+    matParams2.setPe(1e3);
+    matParams2.setLambda(2.0);
+
+    // set liquid saturation
+    fsRef.setSaturation(lPhaseIdx, 0.5);
+    
+    // set pressure of the liquid phase
+    fsRef.setPressure(lPhaseIdx, 1e6);
+
+    // set the remaining parameters of the reference fluid state
+    completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams2, lPhaseIdx);
+  
+    // check the flash calculation
+    checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams2);
+
+    return 0;
+}

test/material/pengrobinson/CMakeLists.txt

+# build target for the simple twophase lens problem
+ADD_EXECUTABLE("test_pengrobinson" test_pengrobinson.cc)
+TARGET_LINK_LIBRARIES("test_pengrobinson" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags 
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+