diff --git a/configure.ac b/configure.ac
index cbaa65656a47ec21a750d38f4fff2b2c801f1c06..485945c6366f38d1e10a32af4c705c3c3a2124c7 100644
--- a/configure.ac
+++ b/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
])
diff --git a/dumux/material/constraintsolvers/immiscibleflash.hh b/dumux/material/constraintsolvers/immiscibleflash.hh
index 6f4f83a148c14b9fc54e1babede6f5c8dd046ed7..0a821c506d507e689280405b002ddb4e67cebd02 100644
--- a/dumux/material/constraintsolvers/immiscibleflash.hh
+++ b/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
diff --git a/dumux/material/fluidsystems/spe5parametercache.hh b/dumux/material/fluidsystems/spe5parametercache.hh
index 53ab2b92a3083441f98c3d2dbd689dfe771e35cc..add88a9b6ddaf3ebde9f9adc9857ae11f284202d 100644
--- a/dumux/material/fluidsystems/spe5parametercache.hh
+++ b/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 };
diff --git a/test/material/CMakeLists.txt b/test/material/CMakeLists.txt
index bfe874093dffc233ebcd83a7f7555c6e8b266df4..3ec7713510173014cc149f83e0e5f21339664b9b 100644
--- a/test/material/CMakeLists.txt
+++ b/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")
diff --git a/test/material/Makefile.am b/test/material/Makefile.am
index 1e361cae2c1b7a7e4b61f6f74ccee6675dcdadb9..b9eb1146bb74c3eb866c748ec12805887a353656 100644
--- a/test/material/Makefile.am
+++ b/test/material/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = fluidsystems ncpflash pengrobinson tabulation
+SUBDIRS = fluidsystems immiscibleflash ncpflash pengrobinson tabulation
include $(top_srcdir)/am/global-rules
diff --git a/test/material/immiscibleflash/CMakeLists.txt b/test/material/immiscibleflash/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..70cb258dd17cc3a2497d96b5c4fde3985ab90e41
--- /dev/null
+++ b/test/material/immiscibleflash/CMakeLists.txt
@@ -0,0 +1,8 @@
+# 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})
+
diff --git a/test/material/immiscibleflash/Makefile.am b/test/material/immiscibleflash/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..5429543d1feb171b3c829f5a246bdf1545533adb
--- /dev/null
+++ b/test/material/immiscibleflash/Makefile.am
@@ -0,0 +1,5 @@
+check_PROGRAMS = test_immiscibleflash
+
+test_immiscibleflash_SOURCES = test_immiscibleflash.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/material/immiscibleflash/test_immiscibleflash.cc b/test/material/immiscibleflash/test_immiscibleflash.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1931a80f1d883413d9d28f2cd117ed08ae0fe37e
--- /dev/null
+++ b/test/material/immiscibleflash/test_immiscibleflash.cc
@@ -0,0 +1,268 @@
+// -*- 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;
+}
diff --git a/test/material/pengrobinson/CMakeLists.txt b/test/material/pengrobinson/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fb1a38c3ac4e75e4ca6bf88a9ee843e65e65f427
--- /dev/null
+++ b/test/material/pengrobinson/CMakeLists.txt
@@ -0,0 +1,8 @@
+# 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})
+