diff --git a/configure.ac b/configure.ac
index 439cb66f3c69b6f9d3b220ca520f907493a16e18..3bc8c43a891eb626b1335c2cae24cd0180e15be4 100644
--- a/configure.ac
+++ b/configure.ac
@@ -145,6 +145,8 @@ AC_CONFIG_FILES([dumux.pc
test/decoupled/1p/Makefile
test/decoupled/2p/Makefile
test/decoupled/2p2c/Makefile
+ test/material/Makefile
+ test/material/tabulation/Makefile
tutorial/Makefile
])
diff --git a/dumux/material/components/component.hh b/dumux/material/components/component.hh
index 133c04f289c5f9bd15e2b26f8b1e8aad3ae65570..7284ffdcca9ba71d3e0eb856e445f6159bdb3d89 100644
--- a/dumux/material/components/component.hh
+++ b/dumux/material/components/component.hh
@@ -39,6 +39,8 @@ template <class Scalar, class Implementation>
class Component
{
public:
+ static const bool isTabulated = false;
+
/*!
* \brief A default routine for initialization, not needed for components and must not be called.
*
diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh
index 02097ad9ccab9b12b071aafb8f6cbd35d1e6fd33..dbf4b99b4de21aae231e2c4af6bf0880d7dd5795 100644
--- a/dumux/material/components/tabulatedcomponent.hh
+++ b/dumux/material/components/tabulatedcomponent.hh
@@ -54,6 +54,8 @@ template <class Scalar, class RawComponent, bool verbose=true>
class TabulatedComponent
{
public:
+ static const bool isTabulated = true;
+
/*!
* \brief Initialize the tables.
*
@@ -88,8 +90,6 @@ public:
gasEnthalpy_ = new Scalar[nTemp_*nPress_];
liquidEnthalpy_ = new Scalar[nTemp_*nPress_];
- gasInternalEnergy_ = new Scalar[nTemp_*nPress_];
- liquidInternalEnergy_ = new Scalar[nTemp_*nPress_];
gasDensity_ = new Scalar[nTemp_*nPress_];
liquidDensity_ = new Scalar[nTemp_*nPress_];
gasViscosity_ = new Scalar[nTemp_*nPress_];
@@ -119,9 +119,6 @@ public:
try { gasEnthalpy_[i] = RawComponent::gasEnthalpy(temperature, pressure); }
catch (NumericalProblem) { gasEnthalpy_[i] = NaN; };
- try { gasInternalEnergy_[i] = RawComponent::gasInternalEnergy(temperature, pressure); }
- catch (NumericalProblem) { gasInternalEnergy_[i] = NaN; };
-
try { gasDensity_[i] = RawComponent::gasDensity(temperature, pressure); }
catch (NumericalProblem) { gasDensity_[i] = NaN; };
@@ -159,9 +156,6 @@ public:
try { liquidEnthalpy_[i] = RawComponent::liquidEnthalpy(temperature, pressure); }
catch (NumericalProblem) { liquidEnthalpy_[i] = NaN; };
- try { liquidInternalEnergy_[i] = RawComponent::liquidInternalEnergy(temperature, pressure); }
- catch (NumericalProblem) { liquidInternalEnergy_[i] = NaN; };
-
try { liquidDensity_[i] = RawComponent::liquidDensity(temperature, pressure); }
catch (NumericalProblem) { liquidDensity_[i] = NaN; };
@@ -261,7 +255,7 @@ public:
}
/*!
- * \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
+ * \brief Specific enthalpy of the liquid \f$\mathrm{[J/kg]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
@@ -286,35 +280,24 @@ public:
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
{
- Scalar result = interpolateGasTP_(gasInternalEnergy_,
- temperature,
- pressure);
- if (std::isnan(result)) {
- printWarning_("gasInternalEnergy", temperature, pressure);
- return RawComponent::gasInternalEnergy(temperature, pressure);
- }
+ Scalar result =
+ gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure);
return result;
}
/*!
- * \brief Specific enthalpy of the liquid \f$\mathrm{[J/kg]}\f$.
+ * \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
{
- Scalar result = interpolateLiquidTP_(liquidInternalEnergy_,
- temperature,
- pressure);
- if (std::isnan(result)) {
- printWarning_("liquidInternalEnergy", temperature, pressure);
- return RawComponent::liquidInternalEnergy(temperature, pressure);
- }
+ Scalar result =
+ liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
return result;
}
-
/*!
* \brief The pressure of gas in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
*
@@ -471,15 +454,15 @@ private:
return std::numeric_limits<Scalar>::quiet_NaN();
}
- unsigned iT = (unsigned) alphaT;
+ unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
alphaT -= iT;
Scalar alphaP1 = pressLiquidIdx_(p, iT);
Scalar alphaP2 = pressLiquidIdx_(p, iT + 1);
- unsigned iP1 = std::min(nPress_ - 2, (unsigned) alphaP1);
+ unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP1));
+ unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP2));
alphaP1 -= iP1;
- unsigned iP2 = std::min(nPress_ - 2, (unsigned) alphaP2);
alphaP2 -= iP2;
return
@@ -499,15 +482,14 @@ private:
return std::numeric_limits<Scalar>::quiet_NaN();
}
- unsigned iT = (unsigned) alphaT;
+ unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
alphaT -= iT;
Scalar alphaP1 = pressGasIdx_(p, iT);
Scalar alphaP2 = pressGasIdx_(p, iT + 1);
-
- unsigned iP1 = std::min(nPress_ - 2, (unsigned) alphaP1);
+ unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP1));
+ unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP2));
alphaP1 -= iP1;
- unsigned iP2 = std::min(nPress_ - 2, (unsigned) alphaP2);
alphaP2 -= iP2;
return
@@ -591,8 +573,9 @@ private:
// returns the index of an entry in a temperature field
static Scalar pressGasIdx_(Scalar pressure, unsigned tempIdx)
{
- Scalar pressMax = maxGasPressure_(tempIdx);
- return (nPress_ - 1)*(pressure - pressMin_)/(pressMax - pressMin_);
+ Scalar pgMin = minGasPressure_(tempIdx);
+ Scalar pgMax = maxGasPressure_(tempIdx);
+ return (nPress_ - 1)*(pressure - pgMin)/(pgMax - pgMin);
}
// returns the index of an entry in a density field
@@ -674,9 +657,6 @@ private:
static Scalar *gasEnthalpy_;
static Scalar *liquidEnthalpy_;
- static Scalar *gasInternalEnergy_;
- static Scalar *liquidInternalEnergy_;
-
static Scalar *gasDensity_;
static Scalar *liquidDensity_;
@@ -725,10 +705,6 @@ Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::gasEnthalpy_;
template <class Scalar, class RawComponent, bool verbose>
Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::liquidEnthalpy_;
template <class Scalar, class RawComponent, bool verbose>
-Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::gasInternalEnergy_;
-template <class Scalar, class RawComponent, bool verbose>
-Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::liquidInternalEnergy_;
-template <class Scalar, class RawComponent, bool verbose>
Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::gasDensity_;
template <class Scalar, class RawComponent, bool verbose>
Scalar* TabulatedComponent<Scalar, RawComponent, verbose>::liquidDensity_;
diff --git a/test/Makefile.am b/test/Makefile.am
index 2ac6b4d5041aecf7f9abd452f02b4166bbce8a63..6d1dc86f37c2adfbbab10b4275b15c2263dc6cbf 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = boxmodels common decoupled
+SUBDIRS = boxmodels common decoupled material
include $(top_srcdir)/am/global-rules
diff --git a/test/material/Makefile.am b/test/material/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..e01d18a4d7bd60fd3e048ae7aa48167ab7f4dd84
--- /dev/null
+++ b/test/material/Makefile.am
@@ -0,0 +1,4 @@
+SUBDIRS = tabulation .
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/test/material/tabulation/Makefile.am b/test/material/tabulation/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..1fbfd4cc696ab9bb05f2e08b8fa996eb8e07eb24
--- /dev/null
+++ b/test/material/tabulation/Makefile.am
@@ -0,0 +1,5 @@
+check_PROGRAMS = test_tabulation
+
+test_tabulation_SOURCES = test_tabulation.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/material/tabulation/test_tabulation.cc b/test/material/tabulation/test_tabulation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4ad0efe4d694781aeffe9badc258e26a45dd035c
--- /dev/null
+++ b/test/material/tabulation/test_tabulation.cc
@@ -0,0 +1,97 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 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. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This is a program to test the tabulation class for of
+ * individual components.
+ *
+ * It either prints "success" or "error", it does not do anything
+ * else.
+ */
+#include "config.h"
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+
+template <class Scalar>
+void isSame(Scalar v, Scalar vRef, Scalar tol=5e-4)
+{
+ if (std::abs( (v - vRef)/vRef ) > tol) {
+ std::cout << "\nerror: " << (v - vRef)/vRef << "\n";
+ exit(1);
+ }
+}
+
+int main()
+{
+ typedef double Scalar;
+ typedef Dumux::H2O<Scalar> IapwsH2O;
+ typedef Dumux::TabulatedComponent<Scalar, IapwsH2O> TabulatedH2O;
+
+ Scalar tempMin = 274.15;
+ Scalar tempMax = 622.15;
+ int nTemp = (int) (tempMax - tempMin)*3/2;
+
+ Scalar pMin = 10.00;
+ Scalar pMax = IapwsH2O::vaporPressure(tempMax*1.1);
+ int nPress = 600;
+
+ std::cout << "Creating tabulation with " << nTemp*nPress << " entries per quantity\n";
+ TabulatedH2O::init(tempMin, tempMax, nTemp,
+ pMin, pMax, nPress);
+
+ std::cout << "Checking tabulation\n";
+
+ int m = nTemp*3;
+ int n = nPress*3;
+ for (int i = 0; i < m; ++i) {
+ Scalar T = tempMin + (tempMax - tempMin)*Scalar(i)/m;
+
+ if (i % std::max(1, m/1000) == 0) {
+ std::cout << Scalar(i)/m*100 << "% done \r";
+ std::cout.flush();
+ }
+
+ isSame(TabulatedH2O::vaporPressure(T),
+ IapwsH2O::vaporPressure(T),
+ 1e-3);
+ for (int j = 0; j < n; ++j) {
+ Scalar p = pMin + (pMax - pMin)*Scalar(j)/n;
+ if (p < IapwsH2O::vaporPressure(T) * 1.03) {
+ Scalar tol = 5e-4;
+ if (p > IapwsH2O::vaporPressure(T))
+ tol = 5e-2;
+ isSame(TabulatedH2O::gasEnthalpy(T,p), IapwsH2O::gasEnthalpy(T,p), tol);
+ isSame(TabulatedH2O::gasInternalEnergy(T,p), IapwsH2O::gasInternalEnergy(T,p), tol);
+ isSame(TabulatedH2O::gasDensity(T,p), IapwsH2O::gasDensity(T,p), tol);
+ isSame(TabulatedH2O::gasViscosity(T,p), IapwsH2O::gasViscosity(T,p), tol);
+ }
+
+ if (p > IapwsH2O::vaporPressure(T) / 1.03) {
+ Scalar tol = 5e-4;
+ if (p < IapwsH2O::vaporPressure(T))
+ tol = 5e-2;
+ isSame(TabulatedH2O::liquidEnthalpy(T,p), IapwsH2O::liquidEnthalpy(T,p), tol);
+ isSame(TabulatedH2O::liquidInternalEnergy(T,p), IapwsH2O::liquidInternalEnergy(T,p), tol);
+ isSame(TabulatedH2O::liquidDensity(T,p), IapwsH2O::liquidDensity(T,p), tol);
+ isSame(TabulatedH2O::liquidViscosity(T,p), IapwsH2O::liquidViscosity(T,p), tol);
+ }
+ }
+ }
+ std::cout << "\nsuccess\n";
+ return 0;
+}