diff --git a/CMakeLists.txt b/CMakeLists.txt
index b5d625411ffb8b362035b15b8472a88de5873bf4..a493dd9bb014ab13ceaf4f272b0436245d9f979c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -199,6 +199,7 @@ add_subdirectory( DataHolderLib )
add_subdirectory( FileIO )
add_subdirectory( GeoLib )
add_subdirectory( InSituLib )
+add_subdirectory( MaterialsLib/Adsorption )
add_subdirectory( MathLib )
add_subdirectory( MeshLib )
add_subdirectory( MeshGeoToolsLib )
diff --git a/MaterialsLib/Adsorption/Adsorption.cpp b/MaterialsLib/Adsorption/Adsorption.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25cb25a79e8d2638bf2f7db23200d226d48f4195
--- /dev/null
+++ b/MaterialsLib/Adsorption/Adsorption.cpp
@@ -0,0 +1,185 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <logog/include/logog.hpp>
+#include "Adsorption.h"
+
+namespace {
+ const double k_rate = 6.0e-3; // to be specified
+
+ template <typename T>
+ T square(const T& v)
+ {
+ return v * v;
+ }
+}
+
+namespace Adsorption
+{
+
+//Saturation pressure for water used in Nunez
+double AdsorptionReaction::getEquilibriumVapourPressure(const double T_Ads)
+{
+ // critical T and p
+ const double Tc = 647.3; // K
+ const double pc = 221.2e5; // Pa
+ // dimensionless T
+ const double Tr = T_Ads/Tc;
+ const double theta = 1. - Tr;
+ // empirical constants
+ const double c[] = {-7.69123,-26.08023,-168.17065,64.23285,-118.96462,4.16717,20.97506,1.0e9,6.0};
+ const double K[] = {c[0]*theta + c[1]*pow(theta,2) + c[2]*pow(theta,3) + c[3]*pow(theta,4) + c[4]*pow(theta,5),
+ 1. + c[5]*theta + c[6]*pow(theta,2)};
+
+ const double exponent = K[0]/(K[1]*Tr) - theta/(c[7]*pow(theta,2) + c[8]);
+ return pc * exp(exponent); // in Pa
+}
+
+// Evaporation enthalpy of water from Nunez
+double AdsorptionReaction::getEvaporationEnthalpy(double T_Ads) // in kJ/kg
+{
+ T_Ads -= 273.15;
+ if (T_Ads <= 10.){
+ const double c[] = {2.50052e3,-2.1068,-3.57500e-1,1.905843e-1,-5.11041e-2,7.52511e-3,-6.14313e-4,2.59674e-5,-4.421e-7};
+ double hv = 0.;
+ for (size_t i=0; i< sizeof(c)/sizeof(c[0]);i++)
+ hv += c[i] * pow(T_Ads,i);
+ return hv;
+ } else if (T_Ads <= 300.){
+ const double c[] = {2.50043e3,-2.35209,1.91685e-4,-1.94824e-5,2.89539e-7,-3.51199e-9,2.06926e-11,-6.4067e-14,8.518e-17,1.558e-20,-1.122e-22};
+ double hv = 0.;
+ for (size_t i=0; i< sizeof(c)/sizeof(c[0]);i++)
+ hv += c[i] * pow(T_Ads,i);
+ return hv;
+ } else {
+ const double c[] = {2.99866e3,-3.1837e-3,-1.566964e1,-2.514e-6,2.045933e-2,1.0389e-8};
+ return ((c[0] + c[2]*T_Ads + c[4]*pow(T_Ads,2))/(1. + c[1]*T_Ads + c[3]*pow(T_Ads,2) + c[5]*pow(T_Ads,3)));
+ }
+}
+
+
+// evaluate specific heat capacity of adsorbate follwing Nunez
+double AdsorptionReaction::getSpecificHeatCapacity(const double T_Ads)
+{
+ const double c[] = {4.224,-3.716e-3,9.351e-5,-7.1786e-7,-9.1266e-9,2.69247e-10,-2.773104e-12,1.553177e-14,-4.982795e-17,8.578e-20,-6.12423e-23};
+ double cp = 0.;
+ for (unsigned i=0; i< sizeof(c)/sizeof(c[0]);i++)
+ cp += c[i] * pow(T_Ads,i);
+ return cp; // kJ/(kg*K)
+}
+
+
+double AdsorptionReaction::getMolarFraction(double xm, double M_this, double M_other)
+{
+ return M_other*xm/(M_other*xm + M_this*(1.0-xm));
+}
+
+
+double AdsorptionReaction::getMassFraction(double xn, double M_this, double M_other)
+{
+ return M_this*xn/(M_this*xn + M_other*(1.0-xn));
+}
+
+
+double AdsorptionReaction::dMolarFraction(double xm, double M_this, double M_other)
+{
+ return M_other * M_this
+ / square(M_other * xm + M_this * (1.0 - xm));
+}
+
+
+double AdsorptionReaction::getReactionRate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading) const
+{
+ const double A = getPotential(p_Ads, T_Ads, M_Ads);
+ double C_eq = getAdsorbateDensity(T_Ads) * characteristicCurve(A);
+ if (C_eq < 0.0) C_eq = 0.0;
+
+ return k_rate * (C_eq - loading); // scaled with mass fraction
+ // this the rate in terms of loading!
+}
+
+void AdsorptionReaction::getDReactionRate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double /*loading*/,
+ std::array<double, 3> &dqdr) const
+{
+ const double A = getPotential(p_Ads, T_Ads, M_Ads);
+ const double p_S = getEquilibriumVapourPressure(T_Ads);
+ const double dAdT = GAS_CONST * log(p_S/p_Ads) / (M_Ads*1.e3);
+ const double dAdp = - GAS_CONST * T_Ads / M_Ads / p_Ads;
+
+ const double W = characteristicCurve(A);
+ const double dWdA = dCharacteristicCurve(A);
+
+ const double rho_Ads = getAdsorbateDensity(T_Ads);
+ const double drhodT = - rho_Ads * getAlphaT(T_Ads);
+
+ dqdr = std::array<double, 3>{{
+ rho_Ads*dWdA*dAdp,
+ drhodT*W + rho_Ads*dWdA*dAdT,
+ -k_rate
+ }};
+}
+
+
+// Evaluate adsorbtion potential A
+double AdsorptionReaction::getPotential(const double p_Ads, double T_Ads, const double M_Ads) const
+{
+ double A = GAS_CONST * T_Ads * log(getEquilibriumVapourPressure(T_Ads)/p_Ads) / (M_Ads*1.e3); // in kJ/kg = J/g
+ return A;
+}
+
+
+double AdsorptionReaction::getLoading(const double rho_curr, const double rho_dry)
+{
+ return rho_curr / rho_dry - 1.0;
+}
+
+
+// Calculate sorption entropy
+double AdsorptionReaction::getEntropy(const double T_Ads, const double A) const
+{
+ const double epsilon = 1.0e-8;
+
+ //* // This change will also change simulation results.
+ const double W_p = characteristicCurve(A+epsilon);
+ const double W_m = characteristicCurve(A-epsilon);
+ const double dAdlnW = 2.0*epsilon/(log(W_p/W_m));
+ // */
+
+ if (W_p <= 0.0 || W_m <= 0.0)
+ {
+ ERR("characteristic curve in negative region (W-, W+): %g, %g", W_m, W_p);
+ return 0.0;
+ }
+
+ return dAdlnW * getAlphaT(T_Ads);
+}
+
+
+//Calculate sorption enthalpy
+double AdsorptionReaction::getEnthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const
+{
+ // TODO [CL] consider using A as obtained from current loading (needs inverse CC A(W)) instead of p_Vapour, T_Vapour
+ const double A = getPotential(p_Ads, T_Ads, M_Ads);
+
+ return (getEvaporationEnthalpy(T_Ads) + A - T_Ads * getEntropy(T_Ads,A))*1000.0; //in J/kg
+}
+
+
+double AdsorptionReaction::getEquilibriumLoading(
+ const double p_Ads, const double T_Ads, const double M_Ads) const
+{
+ const double A = getPotential(p_Ads, T_Ads, M_Ads);
+ return getAdsorbateDensity(T_Ads) * characteristicCurve(A);
+}
+
+
+
+} // namespace Ads
diff --git a/MaterialsLib/Adsorption/Adsorption.h b/MaterialsLib/Adsorption/Adsorption.h
new file mode 100644
index 0000000000000000000000000000000000000000..bdf1cb01202e40527b828cb3520d8862f4c97cac
--- /dev/null
+++ b/MaterialsLib/Adsorption/Adsorption.h
@@ -0,0 +1,91 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef ADSORPTION_H
+#define ADSORPTION_H
+
+#include <array>
+#include <cmath>
+
+#include "Reaction.h"
+
+namespace Adsorption
+{
+
+const double GAS_CONST = 8.3144621;
+
+const double M_N2 = 0.028;
+const double M_H2O = 0.018;
+
+class AdsorptionReaction : public Reaction
+{
+public:
+ // TODO [CL] move those three methods to water properties class
+ static double getEvaporationEnthalpy(const double T_Ads);
+ static double getEquilibriumVapourPressure(const double T_Ads);
+ static double getSpecificHeatCapacity(const double T_Ads); // TODO [CL] why unused?
+
+ static double getMolarFraction(double xm, double M_this, double M_other);
+ static double getMassFraction(double xn, double M_this, double M_other);
+ static double dMolarFraction(double xm, double M_this, double M_other);
+
+ static double getLoading(const double rho_curr, const double rho_dry);
+
+ double getEquilibriumLoading(const double p_Ads, const double T_Ads, const double M_Ads)
+ const override;
+
+ virtual double getEnthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const override;
+ virtual double getReactionRate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading) const override;
+ /**
+ * @brief get_d_reaction_rate
+ * @param p_Ads
+ * @param T_ads
+ * @param M_Ads
+ * @param loading
+ * @param dqdr array containing the differentials wrt: p, T, C
+ */
+ virtual void getDReactionRate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading,
+ std::array<double, 3>& dqdr) const;
+
+protected:
+ virtual double getAdsorbateDensity(const double T_Ads) const = 0;
+ virtual double getAlphaT(const double T_Ads) const = 0;
+ virtual double characteristicCurve(const double A) const = 0;
+ virtual double dCharacteristicCurve(const double A) const = 0;
+
+private:
+ double getPotential(const double p_Ads, const double T_Ads, const double M_Ads) const;
+ double getEntropy(const double T_Ads, const double A) const;
+};
+
+
+inline double curvePolyfrac(const double* coeffs, const double x)
+{
+ // TODO use Horner scheme
+ return ( coeffs[0] + coeffs[2] * x + coeffs[4] * pow(x,2) + coeffs[6] * pow(x,3) )
+ / ( 1.0 + coeffs[1] * x + coeffs[3] * pow(x,2) + coeffs[5] * pow(x,3) );
+}
+
+inline double dCurvePolyfrac(const double* coeffs, const double x)
+{
+ const double x2 = x*x;
+ const double x3 = x2*x;
+ const double u = coeffs[0] + coeffs[2] * x + coeffs[4] * x2 + coeffs[6] * x3;
+ const double du = coeffs[2] + 2.0*coeffs[4] * x + 3.0*coeffs[6] * x2;
+ const double v = 1.0 + coeffs[1] * x + coeffs[3] * x2 + coeffs[5] * x3;
+ const double dv = coeffs[1] + 2.0*coeffs[3] * x + 3.0*coeffs[5] * x2;
+
+ return (du*v - u*dv) / v / v;
+}
+
+}
+
+#endif // ADSORPTION_H
diff --git a/MaterialsLib/Adsorption/CMakeLists.txt b/MaterialsLib/Adsorption/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4ff732ff0ac53b07981c3d0de37343e7ce3ee254
--- /dev/null
+++ b/MaterialsLib/Adsorption/CMakeLists.txt
@@ -0,0 +1,4 @@
+file(GLOB Adsorption_SOURCES *.cpp)
+file(GLOB Adsorption_HEADERS *.h)
+
+add_library(MaterialsLibAdsorption ${Adsorption_HEADERS} ${Adsorption_SOURCES} )
diff --git a/MaterialsLib/Adsorption/Density100MPa.cpp b/MaterialsLib/Adsorption/Density100MPa.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..06241422bfa48faf10553fa8ba70e2c507b3dabd
--- /dev/null
+++ b/MaterialsLib/Adsorption/Density100MPa.cpp
@@ -0,0 +1,63 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "Density100MPa.h"
+
+namespace
+{
+
+// NaX_HighP_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:57
+const double c[] = {
+ 0.3490302932983226, /* a0 */
+ -0.0014061345691831226, /* a1 */
+ -0.0007399303393402753, /* a2 */
+ 5.129318840267485e-09, /* a3 */
+ 5.243619689772646e-07, /* a4 */
+ 6.347011955956523e-10, /* a5 */
+ -9.919599580166727e-11 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double Density100MPa::getAdsorbateDensity(const double T_Ads) const
+{
+ return -0.0013*T_Ads*T_Ads + 0.3529*T_Ads + 1049.2;
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double Density100MPa::getAlphaT(const double T_Ads) const
+{
+ const double rho = -0.0013*T_Ads*T_Ads+0.3529*T_Ads+1049.2;
+ const double drhodT = -0.0026*T_Ads + 0.3529;
+
+ return - drhodT / rho;
+}
+
+// Characteristic curve. Return W (A)
+double Density100MPa::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double Density100MPa::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/Density100MPa.h b/MaterialsLib/Adsorption/Density100MPa.h
new file mode 100644
index 0000000000000000000000000000000000000000..2403e0acca390c9ece1a6f3e8e863d5360ddb2bc
--- /dev/null
+++ b/MaterialsLib/Adsorption/Density100MPa.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITY100MPA_H
+#define MATERIALSLIB_ADSORPTION_DENSITY100MPA_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class Density100MPa : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITY100MPA_H
diff --git a/MaterialsLib/Adsorption/DensityConst.cpp b/MaterialsLib/Adsorption/DensityConst.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..39d097f5afba770cdbe4b9f1664af86a79c033a5
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityConst.cpp
@@ -0,0 +1,60 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityConst.h"
+#include "DensityHauer.h"
+
+namespace
+{
+
+// NaX_Constant_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:20:05
+const double c[] = {
+ 0.3824098506898007, /* a0 */
+ -0.001316857559708455, /* a1 */
+ -0.0007935756090263691, /* a2 */
+ -1.1600036977157845e-07, /* a3 */
+ 5.610354459181838e-07, /* a4 */
+ 7.113664938298873e-10, /* a5 */
+ -1.0668790477629686e-10 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityConst::getAdsorbateDensity(const double /*T_Ads*/) const
+{
+ return rhoWaterHauer(150.0+273.15);
+}
+
+double DensityConst::getAlphaT(const double /*T_Ads*/) const
+{
+ return 0.0;
+}
+
+// Characteristic curve. Return W (A)
+double DensityConst::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); //cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityConst::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityConst.h b/MaterialsLib/Adsorption/DensityConst.h
new file mode 100644
index 0000000000000000000000000000000000000000..646b09dcb1d5ae95f943783d444c7b6989c4c6fb
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityConst.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYCONST_H
+#define MATERIALSLIB_ADSORPTION_DENSITYCONST_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityConst : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYCONST_H
diff --git a/MaterialsLib/Adsorption/DensityCook.cpp b/MaterialsLib/Adsorption/DensityCook.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e798cd16a5bb4831f467f8f07a678965809e9646
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityCook.cpp
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityCook.h"
+
+namespace
+{
+
+// NaX_Dean_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:42
+const double c[] = {
+ 0.3632627555646154, /* a0 */
+ -0.0014090624975800715, /* a1 */
+ -0.0007717609035743321, /* a2 */
+ 5.03634836561135e-09, /* a3 */
+ 5.478509959282738e-07, /* a4 */
+ 6.36458510620815e-10, /* a5 */
+ -1.037977321231462e-10 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityCook::getAdsorbateDensity(const double T_Ads) const
+{
+ return rhoWaterDean(T_Ads);
+}
+
+double DensityCook::getAlphaT(const double T_Ads) const
+{
+ return alphaTWaterDean(T_Ads);
+}
+
+// Characteristic curve. Return W (A)
+double DensityCook::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); //cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; //m^3/kg
+}
+
+double DensityCook::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityCook.h b/MaterialsLib/Adsorption/DensityCook.h
new file mode 100644
index 0000000000000000000000000000000000000000..59d7062fe980fc19223d233454a002777bfd6d58
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityCook.h
@@ -0,0 +1,57 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYCOOK_H
+#define MATERIALSLIB_ADSORPTION_DENSITYCOOK_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityCook : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+inline double rhoWaterDean(const double T_Ads)
+{
+ const double Tcel = T_Ads - 273.15;
+ const double b[] = { 999.9,2.03E-02,-6.16E-03,2.26E-05,-4.68E-08 };
+ if (Tcel <= 100.) {
+ return b[0] + Tcel * (b[1] + Tcel * (b[2] + Tcel * (b[3] + Tcel * b[4]) ) );
+ }
+ else {
+ const double rho_100 = b[0] + b[1]*1.e2 + b[2]*1.e4 + b[3]*1.e6 + b[4]*1.e8;
+ const double aT_100 = -1./rho_100 * (b[1] + 2.*b[2]*1.e2 + 3.*b[3]*1.e4 + 4.*b[4]*1.e6);
+ return rho_100 * (1. - aT_100*(Tcel-100.));
+ }
+}
+
+inline double alphaTWaterDean(const double T_Ads)
+{
+ const double Tcel = T_Ads - 273.15;
+ const double b[] = { 999.9,2.03E-02,-6.16E-03,2.26E-05,-4.68E-08 };
+ if (Tcel <= 100.) {
+ const double r = b[0] + Tcel * (b[1] + Tcel * (b[2] + Tcel * (b[3] + Tcel * b[4]) ) );
+ return -1.0/r * ( b[1] + Tcel * (2.0*b[2] + Tcel * (3.0*b[3] + Tcel * 4.0*b[4]) ) );
+ }
+ else {
+ const double rho_100 = b[0] + b[1]*1.e2 + b[2]*1.e4 + b[3]*1.e6 + b[4]*1.e8;
+ const double aT_100 = -1./rho_100 * (b[1] + 2.*b[2]*1.e2 + 3.*b[3]*1.e4 + 4.*b[4]*1.e6);
+ return aT_100 / (1. - aT_100*(Tcel-100.));
+ }
+}
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYCOOK_H
diff --git a/MaterialsLib/Adsorption/DensityDubinin.cpp b/MaterialsLib/Adsorption/DensityDubinin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8cb1c740306cdbd181414244dde2913c87d95fed
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityDubinin.cpp
@@ -0,0 +1,95 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityDubinin.h"
+#include "DensityCook.h"
+#include "Adsorption.h"
+
+namespace
+{
+
+// NaX_Dubinin_polyfrac_CC.pickle
+// date extracted 2015-06-23 16:47:50 file mtime 2015-06-23 16:47:23
+const double c[] = {
+ 0.3635538371322433, /* a0 */
+ -0.0014521033261199435, /* a1 */
+ -0.0007855160157616825, /* a2 */
+ 4.385666000850253e-08, /* a3 */
+ 5.567776459188524e-07, /* a4 */
+ 6.026002134230559e-10, /* a5 */
+ -1.0477401124006098e-10 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityDubinin::getAdsorbateDensity(const double T_Ads) const
+{
+ const double Tb = 373.1;
+
+ if (T_Ads < Tb) {
+ return rhoWaterDean(T_Ads);
+ } else {
+ const double Tc = 647.3; // K
+ const double pc = 221.2e5; // Pa
+ // boiling point density
+ const double rhob = rhoWaterDean(Tb);
+ // state values
+ const double R = GAS_CONST;
+ const double M = M_H2O;
+ const double b = R * Tc/(8. * pc); // m^3/mol
+ const double rhom = M/b; // kg/m^3
+ const double rho = rhob - (rhob-rhom)/(Tc-Tb)*(T_Ads-Tb);
+ return rho;
+ }
+}
+
+//Thermal expansivity model for water found in the works of Hauer
+double DensityDubinin::getAlphaT(const double T_Ads) const
+{
+ const double Tb = 373.1;
+ if (T_Ads <= Tb) {
+ return alphaTWaterDean(T_Ads);
+ } else {
+ // critical T and p
+ const double Tc = 647.3; // K
+ // const double rhoc = 322.; // kg/m^3
+ const double pc = 221.2e5; // Pa
+ // boiling point density
+ const double rhob = rhoWaterDean(Tb);
+ // state values
+ const double R = GAS_CONST;
+ const double M = M_H2O;
+ const double b = R * Tc/(8. * pc); // m^3/mol
+ const double rhom = M/(b); // kg/m^3
+ const double rho = rhob - (rhob-rhom)/(Tc-Tb)*(T_Ads-Tb);
+ return ((rhob-rhom)/(Tc-Tb)*1./rho);
+ }
+}
+
+// Characteristic curve. Return W (A)
+double DensityDubinin::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityDubinin::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityDubinin.h b/MaterialsLib/Adsorption/DensityDubinin.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ea23344ade4836bc054366bd727364105d28240
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityDubinin.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H
+#define MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityDubinin : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYDUBININ_H
diff --git a/MaterialsLib/Adsorption/DensityHauer.cpp b/MaterialsLib/Adsorption/DensityHauer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..331265b970eb2455f03ad2d965e5f6b756dff60b
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityHauer.cpp
@@ -0,0 +1,63 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityHauer.h"
+
+namespace
+{
+
+// NaX_Hauer_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:19
+const double c[] = {
+ 0.36490158988356747, /* a0 */
+ -0.0013723270478333963, /* a1 */
+ -0.0007655780628099964, /* a2 */
+ -3.353324854315774e-08, /* a3 */
+ 5.424357157710913e-07, /* a4 */
+ 6.613430586648678e-10, /* a5 */
+ -1.0300151379421499e-10 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityHauer::getAdsorbateDensity(const double T_Ads) const
+{
+ return rhoWaterHauer(T_Ads);
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityHauer::getAlphaT(const double T_Ads) const
+{
+ // data like in python script
+ const double T0 = 283.15, alpha0 = 3.781e-4; //K; 1/K
+
+ return alpha0/(1. - alpha0 * (T_Ads-T0)); //in 1/K
+}
+
+// Characteristic curve. Return W (A)
+double DensityHauer::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityHauer::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityHauer.h b/MaterialsLib/Adsorption/DensityHauer.h
new file mode 100644
index 0000000000000000000000000000000000000000..a9f606ad5d1b659b8f2311451287ee0c4647b318
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityHauer.h
@@ -0,0 +1,37 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYHAUER_H
+#define MATERIALSLIB_ADSORPTION_DENSITYHAUER_H
+
+#include "Adsorption.h"
+#include "DensityCook.h"
+
+namespace Adsorption
+{
+
+class DensityHauer : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+inline double rhoWaterHauer(const double T_Ads)
+{
+ // data like in python script
+ const double T0 = 283.15, rho0 = rhoWaterDean(T0), alpha0 = 3.781e-4; // K; kg/m^3; 1/K
+
+ return rho0 * (1. - alpha0 * (T_Ads-T0)); // in kg/m^3
+}
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYHAUER_H
diff --git a/MaterialsLib/Adsorption/DensityLegacy.cpp b/MaterialsLib/Adsorption/DensityLegacy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..407538cb831d8cb8d18a25d7be069502a86e38c6
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityLegacy.cpp
@@ -0,0 +1,64 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityLegacy.h"
+
+namespace
+{
+
+//parameters from least squares fit (experimental data)
+const double c[] = { 0.34102920966608297,
+ -0.0013106032830951296,
+ -0.00060754147575378876,
+ 3.7843404172683339e-07,
+ 4.0107503869519016e-07,
+ 3.1274595098338057e-10,
+ -7.610441241719489e-11
+ };
+
+}
+
+namespace Adsorption
+{
+
+double DensityLegacy::getAdsorbateDensity(const double T_Ads) const
+{
+ //set reference state for adsorbate EOS in Hauer
+ const double T0 = 293.15, rho0 = 998.084, alpha0 = 2.06508e-4; // K; kg/m^3; 1/K
+
+ return (rho0 * (1. - alpha0 * (T_Ads-T0))); // in kg/m^3
+}
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityLegacy::getAlphaT(const double T_Ads) const
+{
+ //set reference state for adsorbate EOS in Hauer
+ const double T0 = 293.15, alpha0 = 2.06508e-4; // K; 1/K
+
+ return (alpha0/(1. - alpha0 * (T_Ads-T0))); // in 1/K
+}
+
+// Characteristic curve. Return W (A)
+double DensityLegacy::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityLegacy::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityLegacy.h b/MaterialsLib/Adsorption/DensityLegacy.h
new file mode 100644
index 0000000000000000000000000000000000000000..053cc050d5e34ffdc7bf0c709e8697901cafd826
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityLegacy.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H
+#define MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityLegacy : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYLEGACY_H
diff --git a/MaterialsLib/Adsorption/DensityMette.cpp b/MaterialsLib/Adsorption/DensityMette.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b997e1355f38137a6f838c2180c9d1259555ec2
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityMette.cpp
@@ -0,0 +1,66 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityMette.h"
+#include "DensityCook.h"
+
+namespace
+{
+// NaX_Mette_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:26
+const double c[] = {
+ 0.36340572890087813, /* a0 */
+ -0.0013449597038375108, /* a1 */
+ -0.0007581210111121073, /* a2 */
+ -7.331279615575401e-08, /* a3 */
+ 5.365656973806218e-07, /* a4 */
+ 6.854673678427112e-10, /* a5 */
+ -1.0197050219481966e-10 /* a6 */
+};
+}
+
+namespace Adsorption
+{
+
+double DensityMette::getAdsorbateDensity(const double T_Ads) const
+{
+ const double T0 = 293.15;
+ const double rho0 = rhoWaterDean(T0);
+ const double alpha20 = alphaTWaterDean(T0);
+ return rho0 / (1. + alpha20*(T_Ads-T0));
+}
+
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityMette::getAlphaT(const double T_Ads) const
+{
+ const double T0 = 293.15;
+ const double alpha20 = alphaTWaterDean(T0);
+ return alpha20 / (1. + alpha20 * (T_Ads-T0));
+}
+
+
+// Characteristic curve. Return W (A)
+double DensityMette::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityMette::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityMette.h b/MaterialsLib/Adsorption/DensityMette.h
new file mode 100644
index 0000000000000000000000000000000000000000..7019fce7d1bcaa324733c6cdd9c2849180545f75
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityMette.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYMETTE_H
+#define MATERIALSLIB_ADSORPTION_DENSITYMETTE_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityMette : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYMETTE_H
diff --git a/MaterialsLib/Adsorption/DensityNunez.cpp b/MaterialsLib/Adsorption/DensityNunez.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ccb7de2b224c55ffcb9404e7f169b55888dbc63
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityNunez.cpp
@@ -0,0 +1,74 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "DensityNunez.h"
+
+namespace
+{
+
+// NaX_Nunez_polyfrac_CC.pickle
+// date extracted 2015-06-23 15:38:35 file mtime 2015-06-23 15:19:34
+const double c[] = {
+ 0.3631900485031771, /* a0 */
+ -0.0014242280940080726, /* a1 */
+ -0.0007751726942386291, /* a2 */
+ 2.1775655036811842e-08, /* a3 */
+ 5.488166913667265e-07, /* a4 */
+ 6.204064716725214e-10, /* a5 */
+ -1.0345385018952998e-10 /* a6 */
+};
+
+}
+
+namespace Adsorption
+{
+
+double DensityNunez::getAdsorbateDensity(const double T_Ads) const
+{
+ // TODO admissable T range: 273.16 K <= T_Ads <= 633.15 K
+ const double a[] = { 1.0644e3,-8.01905,1.445348e-2,-4.19589e-6,-4.5294e-9 };
+ const double b[] = { -8.039e-3,1.8698e-5,-2.3015e-8,2.3809e-11,-1.388e-14 };
+ const double u = a[0] + T_Ads * (a[1] + T_Ads * (a[2] + T_Ads * (a[3] + T_Ads * a[4]) ) );
+ const double v = 1.0 + T_Ads * (b[0] + T_Ads * (b[1] + T_Ads * (b[2] + T_Ads * (b[3] + T_Ads * b[4]) ) ) );
+ return u/v;
+}
+
+
+// Thermal expansivity model for water found in the works of Hauer
+double DensityNunez::getAlphaT(const double T_Ads) const
+{
+ // TODO admissable T range: 273.16 K <= T_Ads <= 633.15 K
+ const double a[] = { 1.0644e3,-8.01905,1.445348e-2,-4.19589e-6,-4.5294e-9 };
+ const double b[] = { -8.039e-3,1.8698e-5,-2.3015e-8,2.3809e-11,-1.388e-14 };
+ const double u = a[0] + T_Ads * (a[1] + T_Ads * (a[2] + T_Ads * (a[3] + T_Ads * a[4]) ) );
+ const double v = 1.0 + T_Ads * (b[0] + T_Ads * (b[1] + T_Ads * (b[2] + T_Ads * (b[3] + T_Ads * b[4]) ) ) );
+ const double du = a[1] + T_Ads * (2.0*a[2] + T_Ads * (3.0*a[3] + T_Ads * 4.0*a[4]) );
+ const double dv = b[0] + T_Ads * (2.0*b[1] + T_Ads * (3.0*b[2] + T_Ads * (4.0*b[3] + T_Ads * 5.0*b[4]) ) );
+ return dv/v - du/u;
+}
+
+
+// Characteristic curve. Return W (A)
+double DensityNunez::characteristicCurve(const double A) const
+{
+ double W = curvePolyfrac(c, A); // cm^3/g
+
+ if (W < 0.0) {
+ W = 0.0; // TODO [CL] debug output
+ }
+
+ return W/1.e3; // m^3/kg
+}
+
+double DensityNunez::dCharacteristicCurve(const double A) const
+{
+ return dCurvePolyfrac(c, A);
+}
+
+}
diff --git a/MaterialsLib/Adsorption/DensityNunez.h b/MaterialsLib/Adsorption/DensityNunez.h
new file mode 100644
index 0000000000000000000000000000000000000000..74fa6702bfaf43fe81f41d291166ba7982755eac
--- /dev/null
+++ b/MaterialsLib/Adsorption/DensityNunez.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H
+#define MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H
+
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+class DensityNunez : public AdsorptionReaction
+{
+public:
+ double getAdsorbateDensity(const double T_Ads) const;
+ double getAlphaT(const double T_Ads) const;
+ double characteristicCurve(const double A) const;
+ double dCharacteristicCurve(const double A) const;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_DENSITYNUNEZ_H
diff --git a/MaterialsLib/Adsorption/Reaction.cpp b/MaterialsLib/Adsorption/Reaction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f08aab83a18c8d3d0ff5ad3984c4614447c7a5e
--- /dev/null
+++ b/MaterialsLib/Adsorption/Reaction.cpp
@@ -0,0 +1,66 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <logog/include/logog.hpp>
+
+#include "Reaction.h"
+
+#include "Density100MPa.h"
+#include "DensityConst.h"
+#include "DensityCook.h"
+#include "DensityDubinin.h"
+#include "DensityHauer.h"
+#include "DensityLegacy.h"
+#include "DensityMette.h"
+#include "DensityNunez.h"
+
+#include "ReactionCaOH2.h"
+#include "ReactionInert.h"
+#include "ReactionSinusoidal.h"
+
+
+namespace Adsorption
+{
+
+std::unique_ptr<Reaction>
+Reaction::
+newInstance(BaseLib::ConfigTree const& conf)
+{
+ auto const type = conf.getConfParam<std::string>("type");
+
+ if (type == "Z13XBF")
+ return std::unique_ptr<Reaction>(new DensityLegacy);
+ else if (type == "Z13XBF_100MPa")
+ return std::unique_ptr<Reaction>(new Density100MPa);
+ else if (type == "Z13XBF_Const")
+ return std::unique_ptr<Reaction>(new DensityConst);
+ else if (type == "Z13XBF_Cook")
+ return std::unique_ptr<Reaction>(new DensityCook);
+ else if (type == "Z13XBF_Dubinin")
+ return std::unique_ptr<Reaction>(new DensityDubinin);
+ else if (type == "Z13XBF_Hauer")
+ return std::unique_ptr<Reaction>(new DensityHauer);
+ else if (type == "Z13XBF_Mette")
+ return std::unique_ptr<Reaction>(new DensityMette);
+ else if (type == "Z13XBF_Nunez")
+ return std::unique_ptr<Reaction>(new DensityNunez);
+ else if (type == "Inert")
+ return std::unique_ptr<Reaction>(new ReactionInert);
+ else if (type == "Sinusoidal")
+ return std::unique_ptr<Reaction>(new ReactionSinusoidal(conf));
+ else if (type == "CaOH2")
+ return std::unique_ptr<Reaction>(new ReactionCaOH2(conf));
+
+ ERR("Unknown reactive system: %s.", type.c_str());
+ std::abort();
+
+ return nullptr;
+}
+
+} // namespace Adsorption
diff --git a/MaterialsLib/Adsorption/Reaction.h b/MaterialsLib/Adsorption/Reaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..b0a882503d4218a88ed3774292bf4d8e381e303e
--- /dev/null
+++ b/MaterialsLib/Adsorption/Reaction.h
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_REACTION_H
+#define MATERIALSLIB_ADSORPTION_REACTION_H
+
+#include <memory>
+
+namespace BaseLib { class ConfigTree; }
+
+namespace Adsorption
+{
+
+class Reaction
+{
+public:
+ static std::unique_ptr<Reaction> newInstance(BaseLib::ConfigTree const& rsys);
+
+ virtual double getEnthalpy(const double p_Ads, const double T_Ads, const double M_Ads) const = 0;
+ virtual double getReactionRate(const double p_Ads, const double T_Ads,
+ const double M_Ads, const double loading) const = 0;
+
+ // TODO get rid of
+ virtual double getEquilibriumLoading(const double, const double, const double) const {
+ return -1.0;
+ }
+
+ virtual ~Reaction() = default;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_REACTION_H
diff --git a/MaterialsLib/Adsorption/ReactionCaOH2.cpp b/MaterialsLib/Adsorption/ReactionCaOH2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1b59cb54defa66ea6baa50be4144f25eb09b3155
--- /dev/null
+++ b/MaterialsLib/Adsorption/ReactionCaOH2.cpp
@@ -0,0 +1,132 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <cassert>
+#include <logog/include/logog.hpp>
+#include "ReactionCaOH2.h"
+#include "Adsorption.h"
+
+namespace Adsorption
+{
+
+const double ReactionCaOH2::_R = GAS_CONST;
+const double ReactionCaOH2::_reaction_enthalpy = -1.12e+05;
+const double ReactionCaOH2::_reaction_entropy = -143.5;
+const double ReactionCaOH2::_M_carrier = M_N2;
+const double ReactionCaOH2::_M_react = M_H2O;
+
+const double ReactionCaOH2::_tol_l = 1e-4;
+const double ReactionCaOH2::_tol_u = 1.0 - 1e-4;
+const double ReactionCaOH2::_tol_rho = 0.1;
+
+const double ReactionCaOH2::rho_low = 1656.0;
+const double ReactionCaOH2::rho_up = 2200.0;
+
+
+double
+ReactionCaOH2::getEnthalpy(const double, const double, const double) const
+{
+ return - _reaction_enthalpy/_M_react;
+}
+
+double
+ReactionCaOH2::getReactionRate(const double, const double, const double, const double) const
+{
+ ERR("get_reaction_rate do not call directly");
+ std::abort();
+ return -1.0;
+}
+
+
+double ReactionCaOH2::getReactionRate(double const solid_density)
+{
+ _rho_s = solid_density;
+ calculateQR();
+ return _qR;
+}
+
+void ReactionCaOH2::updateParam(
+ double T_solid,
+ double p_gas,
+ double x_react,
+ double rho_s_initial)
+{
+ _T_s = T_solid;
+ _p_gas = p_gas / 1e5; // convert Pa to bar
+ _x_react = x_react;
+ _rho_s = rho_s_initial;
+}
+
+void ReactionCaOH2::calculateQR()
+{
+ // Convert mass fraction into mole fraction
+ const double mol_frac_react = AdsorptionReaction::getMolarFraction(_x_react, _M_react, _M_carrier);
+
+ _p_r_g = std::max(mol_frac_react * _p_gas, 1.0e-3); // avoid illdefined log
+ setChemicalEquilibrium();
+ const double dXdt = CaHydration();
+ _qR = (rho_up - rho_low) * dXdt;
+}
+
+// determine equilibrium temperature and pressure according to van't Hoff
+void ReactionCaOH2::setChemicalEquilibrium()
+{
+ _X_D = (_rho_s - rho_up - _tol_rho)/(rho_low - rho_up - 2.0*_tol_rho) ;
+ _X_D = (_X_D < 0.5) ? std::max(_tol_l,_X_D) : std::min(_X_D,_tol_u); // constrain to interval [tol_l;tol_u]
+
+ _X_H = 1.0 - _X_D;
+
+ // calculate equilibrium
+ // using the p_eq to calculate the T_eq - Clausius-Clapeyron
+ _T_eq = (_reaction_enthalpy/_R) / ((_reaction_entropy/_R) + std::log(_p_r_g)); // unit of p in bar
+ // Alternative: Use T_s as T_eq and calculate p_eq - for Schaube kinetics
+ _p_eq = std::exp((_reaction_enthalpy/_R)/_T_s - (_reaction_entropy/_R));
+}
+
+
+double ReactionCaOH2::CaHydration()
+{
+ double dXdt;
+ // step 3, calculate dX/dt
+#ifdef SIMPLE_KINETICS
+ if ( T_s < T_eq ) // hydration - simple model
+#else
+ if ( _p_r_g > _p_eq ) // hydration - Schaube model
+#endif
+ {
+ //X_H = max(tol_l,X_H); //lower tolerance to avoid oscillations at onset of hydration reaction. Set here so that no residual reaction rate occurs at end of hydration.
+#ifdef SIMPLE_KINETICS // this is from P. Schmidt
+ dXdt = -1.0*(1.0-X_H) * (T_s - T_eq) / T_eq * 0.2 * conversion_rate::x_react;
+#else //this is from Schaube
+ if (_X_H == _tol_u || _rho_s == rho_up)
+ dXdt = 0.0;
+ else if ( (_T_eq-_T_s) >= 50.0)
+ dXdt = 13945.0 * exp(-89486.0/_R/_T_s) * std::pow(_p_r_g/_p_eq - 1.0,0.83) * 3.0 * (_X_D) * std::pow(-1.0*log(_X_D),0.666);
+ else
+ dXdt = 1.0004e-34 * exp(5.3332e4/_T_s) * std::pow(_p_r_g, 6.0) * (_X_D);
+#endif
+ }
+ else // dehydration
+ {
+ //X_D = max(tol_l,X_D); //lower tolerance to avoid oscillations at onset of dehydration reaction. Set here so that no residual reaction rate occurs at end of dehydration.
+#ifdef SIMPLE_KINETICS // this is from P. Schmidt
+ dXdt = -1.0* (1.0-X_D) * (T_s - T_eq) / T_eq * 0.05;
+#else
+ if (_X_D == _tol_u || _rho_s == rho_low)
+ dXdt = 0.0;
+ else if (_X_D < 0.2)
+ dXdt = -1.9425e12 * exp( -1.8788e5/_R/_T_s ) * std::pow(1.0-_p_r_g/_p_eq,3.0)*(_X_H);
+ else
+ dXdt = -8.9588e9 * exp( -1.6262e5/_R/_T_s ) * std::pow(1.0-_p_r_g/_p_eq,3.0)*2.0 * std::pow(_X_H, 0.5);
+#endif
+ }
+ return dXdt;
+}
+
+}
diff --git a/MaterialsLib/Adsorption/ReactionCaOH2.h b/MaterialsLib/Adsorption/ReactionCaOH2.h
new file mode 100644
index 0000000000000000000000000000000000000000..d75416d01f3765e557b69b8837c5f6ef5d5c5ef9
--- /dev/null
+++ b/MaterialsLib/Adsorption/ReactionCaOH2.h
@@ -0,0 +1,87 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_REACTIONCAOH2_H
+#define MATERIALSLIB_ADSORPTION_REACTIONCAOH2_H
+
+#include "BaseLib/ConfigTree.h"
+#include "Reaction.h"
+#include "Adsorption.h"
+
+namespace ProcessLib
+{
+template<typename>
+class TESFEMReactionAdaptorCaOH2;
+}
+
+namespace Adsorption
+{
+
+class ReactionCaOH2 final : public Reaction
+{
+public:
+ explicit ReactionCaOH2(BaseLib::ConfigTree const& conf)
+ : _ode_solver_config{conf.getConfSubtree("ode_solver_config")}
+ {}
+
+ double getEnthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+ const double /*M_Ads*/) const override;
+
+ double getReactionRate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
+ const double /*loading*/) const override;
+
+ const BaseLib::ConfigTree& getOdeSolverConfig() const { return _ode_solver_config; }
+
+ // TODO merge with getReactionRate() above
+ double getReactionRate(double const solid_density);
+
+ void updateParam(double T_solid,
+ double _p_gas,
+ double _x_react,
+ double rho_s_initial);
+
+private:
+ void calculateQR();
+ void setChemicalEquilibrium();
+ double CaHydration();
+
+ static const double _R; //!< universal gas constant in J/mol/K
+ double _rho_s; //!< solid phase density
+ double _p_gas; //!< gas phase pressure in unit bar
+ double _p_r_g; //!< pressure of H2O on gas phase
+ double _p_eq = 1.0; //!< equilibrium pressure in bar
+ double _T_eq; //!< equilibrium temperature
+ double _T_s; //!< solid phase temperature
+ double _qR; //!< rate of solid density change
+ double _x_react; //!< mass fraction of water in gas phase
+ double _X_D; //!< mass fraction of dehydration (CaO) in the solid phase
+ double _X_H; //!< mass fraction of hydration in the solid phase
+
+ //! reaction enthalpy in J/mol; negative for exothermic composition reaction
+ static const double _reaction_enthalpy;
+ static const double _reaction_entropy; //!< reaction entropy in J/mol/K
+ static const double _M_carrier; //!< inert component molar mass
+ static const double _M_react; //!< reactive component molar mass
+
+ static const double _tol_l;
+ static const double _tol_u;
+ static const double _tol_rho;
+
+ const BaseLib::ConfigTree _ode_solver_config;
+
+ template<typename>
+ friend class ProcessLib::TESFEMReactionAdaptorCaOH2;
+
+public:
+ static const double rho_low; //! lower density limit
+ static const double rho_up; //! upper density limit
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_REACTIONCAOH2_H
diff --git a/MaterialsLib/Adsorption/ReactionInert.h b/MaterialsLib/Adsorption/ReactionInert.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d53a28eec83381ada00e109182c709c8c673b15
--- /dev/null
+++ b/MaterialsLib/Adsorption/ReactionInert.h
@@ -0,0 +1,37 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_REACTIONINERT_H
+#define MATERIALSLIB_ADSORPTION_REACTIONINERT_H
+
+#include "Reaction.h"
+
+namespace Adsorption
+{
+
+class ReactionInert final : public Reaction
+{
+public:
+ double getEnthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+ const double /*M_Ads*/) const override
+ {
+ return 0.0;
+ }
+
+ double getReactionRate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
+ const double /*loading*/) const override
+ {
+ ERR("Method getReactionRate() should never be called directly");
+ std::abort();
+ return 0.0;
+ }
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_REACTIONINERT_H
diff --git a/MaterialsLib/Adsorption/ReactionSinusoidal.h b/MaterialsLib/Adsorption/ReactionSinusoidal.h
new file mode 100644
index 0000000000000000000000000000000000000000..9fb4bb98abf30daba41ec6ee659776fce5d29b50
--- /dev/null
+++ b/MaterialsLib/Adsorption/ReactionSinusoidal.h
@@ -0,0 +1,48 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALSLIB_ADSORPTION_REACTIONSINUSOIDAL_H
+#define MATERIALSLIB_ADSORPTION_REACTIONSINUSOIDAL_H
+
+#include <logog/include/logog.hpp>
+
+#include "Reaction.h"
+#include "BaseLib/ConfigTree.h"
+
+namespace Adsorption
+{
+
+class ReactionSinusoidal final : public Reaction
+{
+public:
+ explicit ReactionSinusoidal(BaseLib::ConfigTree const& conf)
+ : _enthalpy(conf.getConfParam<double>("reaction_enthalpy"))
+ {
+ }
+
+ double getEnthalpy(const double /*p_Ads*/, const double /*T_Ads*/,
+ const double /*M_Ads*/) const override
+ {
+ return _enthalpy;
+ }
+
+ double getReactionRate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
+ const double /*loading*/) const override
+ {
+ ERR("Method getReactionRate() should never be called directly");
+ std::abort();
+ return 0.0;
+ }
+
+private:
+ double _enthalpy;
+};
+
+}
+#endif // MATERIALSLIB_ADSORPTION_REACTIONSINUSOIDAL_H
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 3a84b340bf711ef059398f4d0a4e6a4e6f141b7e..e1738ab98407c1a7a8e2d555a6a4d91c39854663 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -10,6 +10,7 @@ add_library(ProcessLib ${SOURCES})
target_link_libraries(ProcessLib
AssemblerLib
+ MaterialsLibAdsorption
MeshGeoToolsLib
NumLib # for shape matrices
${VTK_LIBRARIES}