[PL] added TES process

Applications/ApplicationsLib/ProjectData.cpp

@@ -35,6 +35,7 @@
 #include "UncoupledProcessesTimeLoop.h"
 
 #include "ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h"
+#include "ProcessLib/TES/TESProcess.h"
 
 
 namespace detail
@@ -178,6 +179,13 @@ void ProjectData::buildProcesses()
                     *_mesh_vec[0], *nl_slv, std::move(time_disc),
                     _process_variables, _parameters, pc));
         }
+        else if (type == "TES")
+        {
+            _processes.emplace_back(
+                ProcessLib::TES::createTESProcess<GlobalSetupType>(
+                    *_mesh_vec[0], *nl_slv, std::move(time_disc),
+                    _process_variables, _parameters, pc));
+        }
         else
         {
             ERR("Unknown process type: %s", type.c_str());

ProcessLib/CMakeLists.txt

@@ -6,6 +6,9 @@ APPEND_SOURCE_FILES(SOURCES)
 add_subdirectory(GroundwaterFlow)
 APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
 
+add_subdirectory(TES)
+APPEND_SOURCE_FILES(SOURCES TES)
+
 add_library(ProcessLib ${SOURCES})
 
 target_link_libraries(ProcessLib

ProcessLib/TES/.gitignore

+TESLocalAssemblerInner-impl-incl.h

ProcessLib/TES/CMakeLists.txt

+if(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES)
+    configure_file(
+        "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl-dynamic.h.in"
+        "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl.h" COPYONLY
+        )
+else()
+    configure_file(
+        "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl-fixed.h.in"
+        "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl.h" COPYONLY
+        )
+endif()

ProcessLib/TES/TESAssemblyParams.h

+/**
+ * \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 PROCESS_LIB_TESPROCESS_NOTPL_H_
+#define PROCESS_LIB_TESPROCESS_NOTPL_H_
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+
+#include "MaterialsLib/Adsorption/Reaction.h"
+
+#include "ProcessLib/VariableTransformation.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+const unsigned NODAL_DOF = 3;
+const unsigned COMPONENT_ID_PRESSURE = 0;
+const unsigned COMPONENT_ID_TEMPERATURE = 1;
+const unsigned COMPONENT_ID_MASS_FRACTION = 2;
+
+const double M_N2 = 0.028013;
+const double M_H2O = 0.018016;
+
+struct AssemblyParams
+{
+    Trafo trafo_p{1.0};
+    Trafo trafo_T{1.0};
+    Trafo trafo_x{1.0};
+
+    std::unique_ptr<Adsorption::Reaction> react_sys;
+
+    double fluid_specific_heat_source =
+        std::numeric_limits<double>::quiet_NaN();
+    double cpG = std::numeric_limits<
+        double>::quiet_NaN();  // specific isobaric fluid heat capacity
+
+    Eigen::MatrixXd solid_perm_tensor = Eigen::MatrixXd::Constant(
+        3, 3, std::numeric_limits<double>::quiet_NaN());  // TODO get dimensions
+    double solid_specific_heat_source =
+        std::numeric_limits<double>::quiet_NaN();
+    double solid_heat_cond = std::numeric_limits<double>::quiet_NaN();
+    double cpS = std::numeric_limits<
+        double>::quiet_NaN();  // specific isobaric solid heat capacity
+
+    double tortuosity = std::numeric_limits<double>::quiet_NaN();
+    double diffusion_coefficient_component =
+        std::numeric_limits<double>::quiet_NaN();
+
+    double poro = std::numeric_limits<double>::quiet_NaN();
+
+    double rho_SR_dry = std::numeric_limits<double>::quiet_NaN();
+
+    const double M_inert = M_N2;  // N2
+    const double M_react = M_H2O;
+
+    // TODO unify variable names
+    double initial_solid_density = std::numeric_limits<double>::quiet_NaN();
+
+    double delta_t = std::numeric_limits<double>::quiet_NaN();
+    unsigned iteration_in_current_timestep = 0;
+
+    bool output_element_matrices = false;
+
+    unsigned number_of_try_of_iteration = 0;
+    double current_time = std::numeric_limits<double>::quiet_NaN();
+
+    //! Output global matrix/rhs after first iteration.
+    std::size_t timestep = 0;
+    std::size_t total_iteration = 0;
+};
+
+}  // namespace TES
+
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_TESPROCESS_NOTPL_H_

ProcessLib/TES/TESLocalAssembler-impl.h

+/**
+ * \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 PROCESS_LIB_TES_FEM_IMPL_H_
+#define PROCESS_LIB_TES_FEM_IMPL_H_
+
+#include "MaterialsLib/Adsorption/Adsorption.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "TESLocalAssembler.h"
+#include "TESReactionAdaptor.h"
+
+namespace
+{
+enum class MatOutType
+{
+    OGS5,
+    PYTHON
+};
+
+const MatOutType MATRIX_OUTPUT_FORMAT = MatOutType::PYTHON;
+
+// TODO move to some location in the OGS core.
+template <typename Mat>
+void ogs5OutMat(const Mat& mat)
+{
+    for (unsigned r = 0; r < mat.rows(); ++r)
+    {
+        switch (MATRIX_OUTPUT_FORMAT)
+        {
+            case MatOutType::OGS5:
+                if (r != 0)
+                    std::printf("\n");
+                std::printf("|");
+                break;
+            case MatOutType::PYTHON:
+                if (r != 0)
+                    std::printf(",\n");
+                std::printf("[");
+                break;
+        }
+
+        for (unsigned c = 0; c < mat.cols(); ++c)
+        {
+            switch (MATRIX_OUTPUT_FORMAT)
+            {
+                case MatOutType::OGS5:
+                    std::printf(" %.16e", mat(r, c));
+                    break;
+                case MatOutType::PYTHON:
+                    if (c != 0)
+                        std::printf(",");
+                    std::printf(" %23.16g", mat(r, c));
+                    break;
+            }
+        }
+
+        switch (MATRIX_OUTPUT_FORMAT)
+        {
+            case MatOutType::OGS5:
+                std::printf(" | ");
+                break;
+            case MatOutType::PYTHON:
+                std::printf(" ]");
+                break;
+        }
+    }
+    std::printf("\n");
+}
+
+template <typename Vec>
+void ogs5OutVec(const Vec& vec)
+{
+    for (unsigned r = 0; r < vec.size(); ++r)
+    {
+        switch (MATRIX_OUTPUT_FORMAT)
+        {
+            case MatOutType::OGS5:
+                if (r != 0)
+                    std::printf("\n");
+                std::printf("| %.16e | ", vec[r]);
+                break;
+            case MatOutType::PYTHON:
+                if (r != 0)
+                    std::printf(",\n");
+                std::printf("[ %23.16g ]", vec[r]);
+                break;
+        }
+    }
+    std::printf("\n");
+}
+
+}  // anonymous namespace
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+    GlobalDim>::TESLocalAssembler(MeshLib::Element const& e,
+                                  std::size_t const /*local_matrix_size*/,
+                                  unsigned const integration_order,
+                                  AssemblyParams const& asm_params)
+    : _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                        IntegrationMethod_, GlobalDim>(
+          e, integration_order)),
+      _d(asm_params,
+         // TODO narrowing conversion
+         static_cast<const unsigned>(
+             _shape_matrices.front()
+                 .N.rows()) /* number of integration points */,
+         GlobalDim),
+      _local_M(ShapeFunction::NPOINTS * NODAL_DOF,
+               ShapeFunction::NPOINTS * NODAL_DOF),
+      _local_K(ShapeFunction::NPOINTS * NODAL_DOF,
+               ShapeFunction::NPOINTS * NODAL_DOF),
+      _local_b(ShapeFunction::NPOINTS * NODAL_DOF),
+      _integration_order(integration_order)
+{
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
+                       GlobalVector,
+                       GlobalDim>::assemble(const double /*t*/,
+                                            std::vector<double> const& local_x)
+{
+    _local_M.setZero();
+    _local_K.setZero();
+    _local_b.setZero();
+
+    IntegrationMethod_ integration_method(_integration_order);
+    unsigned const n_integration_points = integration_method.getNPoints();
+
+    _d.preEachAssemble();
+
+    for (std::size_t ip(0); ip < n_integration_points; ip++)
+    {
+        auto const& sm = _shape_matrices[ip];
+        auto const& wp = integration_method.getWeightedPoint(ip);
+        auto const weight = wp.getWeight();
+
+        _d.assembleIntegrationPoint(ip, local_x, sm.N, sm.dNdx, sm.J, sm.detJ,
+                                    weight, _local_M, _local_K, _local_b);
+    }
+
+    if (_d.getAssemblyParameters().output_element_matrices)
+    {
+        std::puts("### Element: ?");
+
+        std::puts("---Velocity of water");
+        for (auto const& vs : _d.getData().velocity)
+        {
+            std::printf("| ");
+            for (auto v : vs)
+            {
+                std::printf("%23.16e ", v);
+            }
+            std::printf("|\n");
+        }
+
+        std::printf("\n---Mass matrix: \n");
+        ogs5OutMat(_local_M);
+        std::printf("\n");
+
+        std::printf("---Laplacian + Advective + Content matrix: \n");
+        ogs5OutMat(_local_K);
+        std::printf("\n");
+
+        std::printf("---RHS: \n");
+        ogs5OutVec(_local_b);
+        std::printf("\n");
+    }
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
+                       GlobalVector, GlobalDim>::
+    addToGlobal(
+        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const
+{
+    M.add(indices, _local_M);
+    K.add(indices, _local_K);
+    b.add(indices.rows, _local_b);
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+    GlobalDim>::getIntegrationPointValues(TESIntPtVariables const var,
+                                          std::vector<double>& cache) const
+{
+    return _d.getIntegrationPointValues(var, cache);
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+bool TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+    GlobalDim>::checkBounds(std::vector<double> const& local_x,
+                            std::vector<double> const& local_x_prev_ts)
+{
+    return _d.getReactionAdaptor().checkBounds(local_x, local_x_prev_ts);
+}
+
+}  // namespace TES
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_TES_FEM_IMPL_H_

ProcessLib/TES/TESLocalAssembler.h

+/**
+ * \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 PROCESS_LIB_TES_FEM_H_
+#define PROCESS_LIB_TES_FEM_H_
+
+#include <memory>
+#include <vector>
+
+#include "TESAssemblyParams.h"
+#include "TESLocalAssemblerInner-fwd.h"
+
+#include "NumLib/Extrapolation/Extrapolator.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalMatrix, typename GlobalVector>
+class TESLocalAssemblerInterface
+    : public NumLib::Extrapolatable<GlobalVector, TESIntPtVariables>
+{
+public:
+    virtual ~TESLocalAssemblerInterface() = default;
+
+    virtual void assemble(double const t,
+                          std::vector<double> const& local_x) = 0;
+
+    virtual void addToGlobal(
+        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const = 0;
+
+    virtual bool checkBounds(std::vector<double> const& local_x,
+                             std::vector<double> const& local_x_prev_ts) = 0;
+};
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+class TESLocalAssembler final
+    : public TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>
+{
+public:
+    using ShapeFunction = ShapeFunction_;
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+    TESLocalAssembler(MeshLib::Element const& e,
+                      std::size_t const local_matrix_size,
+                      unsigned const integration_order,
+                      AssemblyParams const& asm_params);
+
+    void assemble(double const t, std::vector<double> const& local_x) override;
+
+    void addToGlobal(
+        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const override;
+
+    Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N = _shape_matrices[integration_point].N;
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::VectorXd>(N.data(), N.size());
+    }
+
+    bool checkBounds(std::vector<double> const& local_x,
+                     std::vector<double> const& local_x_prev_ts) override;
+
+    std::vector<double> const& getIntegrationPointValues(
+        TESIntPtVariables const var, std::vector<double>& cache) const override;
+
+private:
+    std::vector<ShapeMatrices> _shape_matrices;
+
+    using LAT = LocalAssemblerTraits<ShapeMatricesType, ShapeFunction::NPOINTS,
+                                     NODAL_DOF, GlobalDim>;
+
+    TESLocalAssemblerInner<LAT> _d;
+
+    using NodalMatrixType = typename LAT::LocalMatrix;
+    using NodalVectorType = typename LAT::LocalVector;
+
+    static_assert(
+        std::is_same<NodalMatrixType, typename LAT::LocalMatrix>::value,
+        "local matrix and data traits matrix do not coincide");
+    static_assert(
+        std::is_same<NodalVectorType, typename LAT::LocalVector>::value,
+        "local vector and data traits vector do not coincide");
+
+    // TODO Change VectorMatrixAssembler s.t. these can be omitted.
+    NodalMatrixType _local_M;
+    NodalMatrixType _local_K;
+    NodalVectorType _local_b;
+
+    // TODO Use the value from Process
+    unsigned const _integration_order;
+};
+
+}  // namespace TES
+}  // namespace ProcessLib
+
+#include "TESLocalAssembler-impl.h"
+
+#endif  // PROCESS_LIB_TES_FEM_H_

ProcessLib/TES/TESLocalAssemblerData.cpp

+/**
+ * \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 "TESLocalAssemblerData.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+TESLocalAssemblerData::TESLocalAssemblerData(AssemblyParams const& ap_,
+                                             const unsigned num_int_pts,
+                                             const unsigned dimension)
+    : ap(ap_),
+      solid_density(num_int_pts, ap_.initial_solid_density),
+      reaction_rate(num_int_pts),
+      velocity(dimension, std::vector<double>(num_int_pts)),
+      reaction_adaptor(TESFEMReactionAdaptor::newInstance(*this)),
+      solid_density_prev_ts(num_int_pts, ap_.initial_solid_density),
+      reaction_rate_prev_ts(num_int_pts)
+{
+}
+
+TESLocalAssemblerData::~TESLocalAssemblerData() = default;
+}
+}  // namespaces

ProcessLib/TES/TESLocalAssemblerData.h

+/**
+ * \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 PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H
+#define PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H
+
+#include "TESAssemblyParams.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+class TESFEMReactionAdaptor;
+
+struct TESLocalAssemblerData
+{
+    TESLocalAssemblerData(AssemblyParams const& ap_, const unsigned num_int_pts,
+                          const unsigned dimension);
+
+    ~TESLocalAssemblerData();
+
+    AssemblyParams const& ap;
+
+    // integration point quantities
+    std::vector<double> solid_density;
+    std::vector<double> reaction_rate;  // dC/dt * rho_SR_dry
+    std::vector<std::vector<double>>
+        velocity;  // vector of velocities for each integration point
+
+    // integration point values of unknowns -- temporary storage
+    double p = std::numeric_limits<double>::quiet_NaN();  // gas pressure
+    double T = std::numeric_limits<double>::quiet_NaN();  // temperature
+    double vapour_mass_fraction = std::numeric_limits<
+        double>::quiet_NaN();  // fluid mass fraction of the second component
+
+    // temporary storage for some properties
+    // values do not change during the assembly of one integration point
+    double rho_GR = std::numeric_limits<double>::quiet_NaN();
+    double p_V =
+        std::numeric_limits<double>::quiet_NaN();  // vapour partial pressure
+    double qR = std::numeric_limits<
+        double>::quiet_NaN();  // reaction rate, use this in assembly!!!
+
+    std::unique_ptr<TESFEMReactionAdaptor> const reaction_adaptor;
+
+    // variables at previous timestep
+    std::vector<double> solid_density_prev_ts;
+    std::vector<double> reaction_rate_prev_ts;  // could also be calculated from
+                                                // solid_density_prev_ts
+};
+}
+}  // namespaces
+#endif  // PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H

ProcessLib/TES/TESLocalAssemblerInner-fwd.h

+/**
+ * \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 PROCESS_LIB_TESFEM_DATA_FWD_H_
+#define PROCESS_LIB_TESFEM_DATA_FWD_H_
+
+#include "TESLocalAssemblerInner.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+
+extern template class TESLocalAssemblerInner<
+    LocalAssemblerTraits<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
+              "inconsistent use of macros");
+#else
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
+              "inconsistent use of macros");
+#endif
+}
+}
+
+#endif  // PROCESS_LIB_TESFEM_DATA_FWD_H_

ProcessLib/TES/TESLocalAssemblerInner-impl-incl-dynamic.h.in

+#ifndef PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+#define PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1, "inconsistent use of macros");
+
+#endif // PROCESSLIB_TES_FEM_DATA_IMPL_INCL

ProcessLib/TES/TESLocalAssemblerInner-impl-incl-fixed.h.in

+#ifndef PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+#define PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0, "inconsistent use of macros");
+#include "TESLocalAssemblerInner-impl.h"
+
+#endif // PROCESSLIB_TES_FEM_DATA_IMPL_INCL

ProcessLib/TES/TESLocalAssemblerInner-impl.h

+/**
+ * \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
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#ifndef PROCESS_LIB_TESFEM_DATA_IMPL_H_
+#define PROCESS_LIB_TESFEM_DATA_IMPL_H_
+
+#include <cstdio>
+#include <iostream>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/Function/Interpolation.h"
+
+#include "TESLocalAssemblerInner-fwd.h"
+#include "TESOGS5MaterialModels.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename Traits>
+TESLocalAssemblerInner<Traits>::TESLocalAssemblerInner(
+    const AssemblyParams& ap, const unsigned num_int_pts,
+    const unsigned dimension)
+    : _d{ap, num_int_pts, dimension}
+{
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getMassCoeffMatrix(
+    const unsigned int_pt)
+{
+    // TODO: Dalton's law property
+    const double dxn_dxm = Adsorption::AdsorptionReaction::dMolarFraction(
+        _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
+
+    const double M_pp = _d.ap.poro / _d.p * _d.rho_GR;
+    const double M_pT = -_d.ap.poro / _d.T * _d.rho_GR;
+    const double M_px = (_d.ap.M_react - _d.ap.M_inert) * _d.p /
+                        (GAS_CONST * _d.T) * dxn_dxm * _d.ap.poro;
+
+    const double M_Tp = -_d.ap.poro;
+    const double M_TT =
+        _d.ap.poro * _d.rho_GR * _d.ap.cpG  // TODO: vapour heat capacity
+        +
+        (1.0 - _d.ap.poro) * _d.solid_density[int_pt] *
+            _d.ap.cpS;  // TODO: adsorbate heat capacity
+    const double M_Tx = 0.0;
+
+    const double M_xp = 0.0;
+    const double M_xT = 0.0;
+    const double M_xx = _d.ap.poro * _d.rho_GR;
+
+    Eigen::Matrix3d M;
+    M << M_pp, M_pT, M_px, M_Tp, M_TT, M_Tx, M_xp, M_xT, M_xx;
+
+    return M;
+}
+
+template <typename Traits>
+typename Traits::LaplaceMatrix
+TESLocalAssemblerInner<Traits>::getLaplaceCoeffMatrix(const unsigned /*int_pt*/,
+                                                      const unsigned dim)
+{
+    const double eta_GR = fluid_viscosity(_d.p, _d.T, _d.vapour_mass_fraction);
+
+    const double lambda_F =
+        fluid_heat_conductivity(_d.p, _d.T, _d.vapour_mass_fraction);
+    const double lambda_S = _d.ap.solid_heat_cond;
+
+    using Mat = typename Traits::MatrixDimDim;
+
+    typename Traits::LaplaceMatrix L =
+        Traits::LaplaceMatrix::Zero(dim * NODAL_DOF, dim * NODAL_DOF);
+
+    // TODO: k_rel
+    // L_pp
+    Traits::blockDimDim(L, 0, 0, dim, dim) =
+        Traits::blockDimDim(_d.ap.solid_perm_tensor, 0, 0, dim, dim) *
+        _d.rho_GR / eta_GR;
+
+    // TODO: add zeolite part
+    // L_TT
+    Traits::blockDimDim(L, dim, dim, dim, dim) =
+        Mat::Identity(dim, dim) *
+        (_d.ap.poro * lambda_F + (1.0 - _d.ap.poro) * lambda_S);
+
+    // L_xx
+    Traits::blockDimDim(L, 2 * dim, 2 * dim, dim, dim) =
+        Mat::Identity(dim, dim) * (_d.ap.tortuosity * _d.ap.poro * _d.rho_GR *
+                                   _d.ap.diffusion_coefficient_component);
+
+    return L;
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getAdvectionCoeffMatrix(
+    const unsigned /*int_pt*/)
+{
+    const double A_pp = 0.0;
+    const double A_pT = 0.0;
+
+    const double A_px = 0.0;
+
+    const double A_Tp = 0.0;
+
+    const double A_TT = _d.rho_GR * _d.ap.cpG;  // porosity?
+    const double A_Tx = 0.0;
+
+    const double A_xp = 0.0;
+    const double A_xT = 0.0;
+    const double A_xx = _d.rho_GR;  // porosity?
+
+    Eigen::Matrix3d A;
+    A << A_pp, A_pT, A_px, A_Tp, A_TT, A_Tx, A_xp, A_xT, A_xx;
+
+    return A;
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getContentCoeffMatrix(
+    const unsigned /*int_pt*/)
+{
+    const double C_pp = 0.0;
+    const double C_pT = 0.0;
+
+    const double C_px = 0.0;
+
+    const double C_Tp = 0.0;
+
+    const double C_TT = 0.0;
+    const double C_Tx = 0.0;
+
+    const double C_xp = 0.0;
+    const double C_xT = 0.0;
+    const double C_xx = (_d.ap.poro - 1.0) * _d.qR;
+
+    Eigen::Matrix3d C;
+    C << C_pp, C_pT, C_px, C_Tp, C_TT, C_Tx, C_xp, C_xT, C_xx;
+
+    return C;
+}
+
+template <typename Traits>
+Eigen::Vector3d TESLocalAssemblerInner<Traits>::getRHSCoeffVector(
+    const unsigned int_pt)
+{
+    const double reaction_enthalpy =
+        _d.ap.react_sys->getEnthalpy(_d.p_V, _d.T, _d.ap.M_react);
+
+    const double rhs_p =
+        (_d.ap.poro - 1.0) * _d.qR;  // TODO [CL] body force term
+
+    const double rhs_T =
+        _d.rho_GR * _d.ap.poro * _d.ap.fluid_specific_heat_source +
+        (1.0 - _d.ap.poro) * _d.qR * reaction_enthalpy +
+        _d.solid_density[int_pt] * (1.0 - _d.ap.poro) *
+            _d.ap.solid_specific_heat_source;
+    // TODO [CL] momentum production term
+
+    const double rhs_x =
+        (_d.ap.poro - 1.0) * _d.qR;  // TODO [CL] what if x < 0.0
+
+    Eigen::Vector3d rhs;
+    rhs << rhs_p, rhs_T, rhs_x;
+
+    return rhs;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::initReaction(const unsigned int_pt)
+{
+    auto const& rate = _d.reaction_adaptor->initReaction(int_pt);
+    _d.qR = rate.reaction_rate;
+    _d.reaction_rate[int_pt] = rate.reaction_rate;
+    _d.solid_density[int_pt] = rate.solid_density;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
+    const unsigned int_pt,
+    const std::vector<double>& localX,
+    typename Traits::ShapeMatrices::ShapeType const& smN,
+    typename Traits::ShapeMatrices::DxShapeType const& /*smDNdx*/,
+    typename Traits::ShapeMatrices::JacobianType const& /*smJ*/,
+    const double /*smDetJ*/)
+{
+#ifndef NDEBUG
+    // fill local data with garbage to aid in debugging
+    _d.p = _d.T = _d.vapour_mass_fraction = _d.p_V = _d.rho_GR = _d.qR =
+        std::numeric_limits<double>::quiet_NaN();
+#endif
+
+    NumLib::shapeFunctionInterpolate(localX, smN, _d.p, _d.T,
+                                     _d.vapour_mass_fraction);
+
+    // pre-compute certain properties
+    _d.p_V = _d.p * Adsorption::AdsorptionReaction::getMolarFraction(
+                        _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
+
+    initReaction(int_pt);
+
+    assert(_d.p > 0.0);
+    assert(_d.T > 0.0);
+    assert(0.0 <= _d.vapour_mass_fraction && _d.vapour_mass_fraction <= 1.0);
+
+    _d.rho_GR = fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
+}
+
+template <typename Traits>
+std::vector<double> const&
+TESLocalAssemblerInner<Traits>::getIntegrationPointValues(
+    TESIntPtVariables const var, std::vector<double>& cache) const
+{
+    switch (var)
+    {
+        case TESIntPtVariables::REACTION_RATE:
+            return _d.reaction_rate;
+        case TESIntPtVariables::SOLID_DENSITY:
+            return _d.solid_density;
+        case TESIntPtVariables::VELOCITY_X:
+            return _d.velocity[0];
+        case TESIntPtVariables::VELOCITY_Y:
+            assert(_d.velocity.size() >= 2);
+            return _d.velocity[1];
+        case TESIntPtVariables::VELOCITY_Z:
+            assert(_d.velocity.size() >= 3);
+            return _d.velocity[2];
+
+        case TESIntPtVariables::LOADING:
+        {
+            auto& Cs = cache;
+            Cs.clear();
+            Cs.reserve(_d.solid_density.size());
+
+            for (auto rho_SR : _d.solid_density)
+            {
+                Cs.push_back(rho_SR / _d.ap.rho_SR_dry - 1.0);
+            }
+
+            return Cs;
+        }
+
+        // TODO that's an element value, ain't it?
+        case TESIntPtVariables::REACTION_DAMPING_FACTOR:
+        {
+            auto const num_integration_points = _d.solid_density.size();
+            auto& alphas = cache;
+            alphas.clear();
+            alphas.resize(num_integration_points,
+                          _d.reaction_adaptor->getReactionDampingFactor());
+
+            return alphas;
+        }
+    }
+
+    cache.clear();
+    return cache;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
+    unsigned integration_point,
+    std::vector<double> const& localX,
+    const typename Traits::ShapeMatrices::ShapeType& smN,
+    const typename Traits::ShapeMatrices::DxShapeType& smDNdx,
+    const typename Traits::ShapeMatrices::JacobianType& smJ,
+    const double smDetJ,
+    const double weight,
+    typename Traits::LocalMatrix& local_M,
+    typename Traits::LocalMatrix& local_K,
+    typename Traits::LocalVector& local_b)
+{
+    preEachAssembleIntegrationPoint(integration_point, localX, smN, smDNdx, smJ,
+                                    smDetJ);
+
+    auto const N = smDNdx.cols();  // number of integration points
+    auto const D = smDNdx.rows();  // global dimension: 1, 2 or 3
+
+    assert(N * NODAL_DOF == local_M.cols());
+
+    auto const laplaceCoeffMat = getLaplaceCoeffMatrix(integration_point, D);
+    assert(laplaceCoeffMat.cols() == D * NODAL_DOF);
+    auto const massCoeffMat = getMassCoeffMatrix(integration_point);
+    auto const advCoeffMat = getAdvectionCoeffMatrix(integration_point);
+    auto const contentCoeffMat = getContentCoeffMatrix(integration_point);
+
+    // calculate velocity
+    assert((unsigned)smDNdx.rows() == _d.velocity.size() &&
+           (unsigned)smDNdx.cols() == _d.velocity[0].size());
+
+    auto const velocity =
+        (Traits::blockDimDim(laplaceCoeffMat, 0, 0, D, D) *
+         (smDNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
+                                                                  N)  // grad_p
+          /
+          -_d.rho_GR))
+            .eval();
+    assert(velocity.size() == D);
+
+    for (unsigned d = 0; d < D; ++d)
+    {
+        _d.velocity[d][integration_point] = velocity[d];
+    }
+
+    auto const detJ_w_N = (smDetJ * weight * smN).eval();
+    auto const detJ_w_N_NT = (detJ_w_N * smN.transpose()).eval();
+    assert(detJ_w_N_NT.rows() == N && detJ_w_N_NT.cols() == N);
+
+    auto const detJ_w_N_vT_dNdx =
+        (detJ_w_N * velocity.transpose() * smDNdx).eval();
+    assert(detJ_w_N_vT_dNdx.rows() == N && detJ_w_N_vT_dNdx.cols() == N);
+
+    for (unsigned r = 0; r < NODAL_DOF; ++r)
+    {
+        for (unsigned c = 0; c < NODAL_DOF; ++c)
+        {
+            Traits::blockShpShp(local_K, N * r, N * c, N, N).noalias() +=
+                smDetJ * weight * smDNdx.transpose() *
+                    Traits::blockDimDim(laplaceCoeffMat, D * r, D * c, D, D) *
+                    smDNdx  // end Laplacian part
+                + detJ_w_N_NT * contentCoeffMat(r, c) +
+                detJ_w_N_vT_dNdx * advCoeffMat(r, c);
+            Traits::blockShpShp(local_M, N * r, N * c, N, N).noalias() +=
+                detJ_w_N_NT * massCoeffMat(r, c);
+        }
+    }
+
+    auto const rhsCoeffVector = getRHSCoeffVector(integration_point);
+
+    for (unsigned r = 0; r < NODAL_DOF; ++r)
+    {
+        Traits::blockShp(local_b, N * r, N).noalias() +=
+            rhsCoeffVector(r) * smN * smDetJ * weight;
+    }
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::preEachAssemble()
+{
+    if (_d.ap.iteration_in_current_timestep == 1)
+    {
+        if (_d.ap.number_of_try_of_iteration ==
+            1)  // TODO has to hold if the above holds.
+        {
+            _d.solid_density_prev_ts = _d.solid_density;
+            _d.reaction_rate_prev_ts = _d.reaction_rate;
+
+            _d.reaction_adaptor->preZerothTryAssemble();
+        }
+        else
+        {
+            _d.solid_density = _d.solid_density_prev_ts;
+        }
+    }
+}
+
+}  // namespace TES
+
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_TESFEM_DATA_IMPL_H_

ProcessLib/TES/TESLocalAssemblerInner.cpp

+/**
+ * \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
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#include "TESLocalAssemblerInner-fwd.h"
+
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+#include "TESLocalAssemblerInner-impl.h"
+#endif
+
+namespace ProcessLib
+{
+namespace TES
+{
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+
+template class TESLocalAssemblerInner<
+    LocalAssemblerTraits<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
+              "inconsistent use of macros");
+#else
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
+              "inconsistent use of macros");
+#endif
+}
+}

ProcessLib/TES/TESLocalAssemblerInner.h

+/**
+ * \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
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#ifndef PROCESS_LIB_TES_FEM_NOTPL_H_
+#define PROCESS_LIB_TES_FEM_NOTPL_H_
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/VariableTransformation.h"
+
+#include "TESLocalAssemblerData.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+enum class TESIntPtVariables : unsigned
+{
+    SOLID_DENSITY,
+    REACTION_RATE,
+    VELOCITY_X,
+    VELOCITY_Y,
+    VELOCITY_Z,
+    LOADING,
+    REACTION_DAMPING_FACTOR
+};
+
+template <typename Traits>
+class TESLocalAssemblerInner
+{
+public:
+    explicit TESLocalAssemblerInner(AssemblyParams const& ap,
+                                    const unsigned num_int_pts,
+                                    const unsigned dimension);
+
+    void assembleIntegrationPoint(
+        unsigned integration_point,
+        std::vector<double> const& localX,
+        typename Traits::ShapeMatrices::ShapeType const& smN,
+        typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
+        typename Traits::ShapeMatrices::JacobianType const& smJ,
+        const double smDetJ,
+        const double weight,
+        typename Traits::LocalMatrix& local_M,
+        typename Traits::LocalMatrix& local_K,
+        typename Traits::LocalVector& local_b);
+
+    void preEachAssemble();
+
+    std::vector<double> const& getIntegrationPointValues(
+        TESIntPtVariables const var, std::vector<double>& cache) const;
+
+    // TODO better encapsulation
+    AssemblyParams const& getAssemblyParameters() const { return _d.ap; }
+    TESFEMReactionAdaptor const& getReactionAdaptor() const
+    {
+        return *_d.reaction_adaptor;
+    }
+    TESFEMReactionAdaptor& getReactionAdaptor() { return *_d.reaction_adaptor; }
+    TESLocalAssemblerData const& getData() const { return _d; }
+private:
+    Eigen::Matrix3d getMassCoeffMatrix(const unsigned int_pt);
+    typename Traits::LaplaceMatrix getLaplaceCoeffMatrix(const unsigned int_pt,
+                                                         const unsigned dim);
+    Eigen::Matrix3d getAdvectionCoeffMatrix(const unsigned int_pt);
+    Eigen::Matrix3d getContentCoeffMatrix(const unsigned int_pt);
+    Eigen::Vector3d getRHSCoeffVector(const unsigned int_pt);
+
+    void preEachAssembleIntegrationPoint(
+        const unsigned int_pt,
+        std::vector<double> const& localX,
+        typename Traits::ShapeMatrices::ShapeType const& smN,
+        typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
+        typename Traits::ShapeMatrices::JacobianType const& smJ,
+        const double smDetJ);
+
+    void initReaction(const unsigned int_pt);
+
+    TESLocalAssemblerData _d;
+};
+
+}  // namespace TES
+
+}  // namespace ProcessLib
+
+// tricking cmake dependency checker
+#include "TESLocalAssemblerInner-impl-incl.h"
+
+#endif  // PROCESS_LIB_TES_FEM_NOTPL_H_

ProcessLib/TES/TESOGS5MaterialModels.cpp

+/**
+ * \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 "TESOGS5MaterialModels.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+const double FluidHeatConductivityN2::A[5] = {0.46649, -0.57015, 0.19164,
+                                              -0.03708, 0.00241};
+
+const double FluidHeatConductivityN2::f[9] = {
+    -0.837079888737e3,   0.37914711487e2,    -0.601737844275,
+    0.350418363823e1,    -0.874955653028e-5, 0.148968607239e-7,
+    -0.256370354277e-11, 0.100773735767e1,   0.335340610e4};
+
+const double FluidHeatConductivityN2::C[4] = {3.3373542, 0.37098251, 0.89913456,
+                                              0.16972505};
+
+const double FluidViscosityN2::A[5] = {0.46649, -0.57015, 0.19164, -0.03708,
+                                       0.00241};
+const double FluidViscosityN2::C[5] = {-20.09997, 3.4376416, -1.4470051,
+                                       -0.027766561, -0.21662362};
+
+const double FluidHeatConductivityH2O::a[4] = {0.0102811, 0.0299621, 0.0156146,
+                                               -0.00422464};
+
+}  // TES
+}  // ProcessLib

ProcessLib/TES/TESOGS5MaterialModels.h

+/**
+ * \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 PROCESSLIB_TES_OGS5MATERIALMODELS
+#define PROCESSLIB_TES_OGS5MATERIALMODELS
+
+#include "TESAssemblyParams.h"
+
+namespace
+{
+const double GAS_CONST = 8.3144621;
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+inline double fluid_density(const double p, const double T, const double x)
+{
+    // OGS-5 density model 26
+
+    const double M0 = ProcessLib::TES::M_N2;
+    const double M1 = ProcessLib::TES::M_H2O;
+
+    const double xn = M0 * x / (M0 * x + M1 * (1.0 - x));
+
+    return p / (GAS_CONST * T) * (M1 * xn + M0 * (1.0 - xn));
+    ;
+}
+
+template <int i>
+double mypow(const double x)
+{
+    if (i < 0)
+    {
+        return 1.0 / mypow<-i>(x);
+    }
+    else
+    {
+        const double p = mypow<(i >> 1)>(x);
+        return (i & 1) ? p * p * x : p * p;
+    }
+}
+
+template <>
+inline double mypow<0>(const double /*x*/)
+{
+    return 1.0;
+}
+
+struct FluidViscosityN2
+{
+    static double get(double rho, double T)
+    {
+        const double rho_c = 314;   // [kg/m3]
+        const double CVF = 14.058;  // [1e-3 Pa-s]
+
+        const double sigma = 0.36502496e-09;
+        const double k = 1.38062e-23;
+        const double eps = 138.08483e-23;
+        const double c1 = 0.3125;
+        const double c2 = 2.0442e-49;
+
+        const double T_star = T * k / eps;
+        rho = rho / rho_c;
+
+        double Omega = loop1_term<0>(T_star);
+        Omega += loop1_term<1>(T_star);
+        Omega += loop1_term<2>(T_star);
+        Omega += loop1_term<3>(T_star);
+        Omega += loop1_term<4>(T_star);
+
+        Omega = std::exp(Omega);
+
+        // eta in [Pa*s]
+        const double eta_0 = c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega);
+
+        double sum = loop2_term<2>(rho);
+        sum += loop2_term<3>(rho);
+        sum += loop2_term<4>(rho);
+
+        //
+        const double eta_r =
+            CVF * 1e-6 * (C[0] / (rho - C[1]) + C[0] / C[1] + sum);
+
+        return eta_0 + eta_r;  // [Pa*s]
+    }
+
+private:
+    template <unsigned i>
+    static double loop1_term(double T_star)
+    {
+        return A[i] * mypow<i>(log(T_star));
+    }
+
+    template <unsigned i>
+    static double loop2_term(double rho)
+    {
+        return C[i] * mypow<i - 1>(rho);
+    }
+
+    static const double A[5];
+    static const double C[5];
+};
+
+struct FluidViscosityH2O
+{
+    static double get(double rho, double T)
+    {
+        double my, my_0, my_1;
+        double H[4];
+
+        T = T / 647.096;
+        rho = rho / 322.0;
+
+        H[0] = 1.67752;
+        H[1] = 2.20462;
+        H[2] = 0.6366564;
+        H[3] = -0.241605;
+
+        double h[6][7] = {0.0};
+        h[0][0] = 0.520094000;
+        h[1][0] = 0.085089500;
+        h[2][0] = -1.083740000;
+        h[3][0] = -0.289555000;
+        h[0][1] = 0.222531000;
+        h[1][1] = 0.999115000;
+        h[2][1] = 1.887970000;
+        h[3][1] = 1.266130000;
+        h[5][1] = 0.120573000;
+        h[0][2] = -0.281378000;
+        h[1][2] = -0.906851000;
+        h[2][2] = -0.772479000;
+        h[3][2] = -0.489837000;
+        h[4][2] = -0.257040000;
+        h[0][3] = 0.161913000;
+        h[1][3] = 0.257399000;
+        h[0][4] = -0.032537200;
+        h[3][4] = 0.069845200;
+        h[4][5] = 0.008721020;
+        h[3][6] = -0.004356730;
+        h[5][6] = -0.000593264;
+
+        double sum1 = H[0] / mypow<0>(T);
+        sum1 += H[1] / mypow<1>(T);
+        sum1 += H[2] / mypow<2>(T);
+        sum1 += H[3] / mypow<3>(T);
+
+        my_0 = 100 * std::sqrt(T) / sum1;
+
+        double sum2 = inner_loop<0>(rho, T, h);
+        sum2 += inner_loop<1>(rho, T, h);
+        sum2 += inner_loop<2>(rho, T, h);
+        sum2 += inner_loop<3>(rho, T, h);
+        sum2 += inner_loop<4>(rho, T, h);
+        sum2 += inner_loop<5>(rho, T, h);
+
+        my_1 = std::exp(rho * sum2);
+
+        my = (my_0 * my_1) / 1e6;
+        return my;
+    }
+
+private:
+    template <int i>
+    static double inner_loop(const double rho,
+                             const double T,
+                             const double (&h)[6][7])
+    {
+        const double base = rho - 1.0;
+
+        double sum3 = h[i][0] * mypow<0>(base);
+        sum3 += h[i][1] * mypow<1>(base);
+        sum3 += h[i][2] * mypow<2>(base);
+        sum3 += h[i][3] * mypow<3>(base);
+        sum3 += h[i][4] * mypow<4>(base);
+        sum3 += h[i][5] * mypow<5>(base);
+        sum3 += h[i][6] * mypow<6>(base);
+
+        return mypow<i>(1 / T - 1) * sum3;
+    }
+};
+
+inline double fluid_viscosity(const double p, const double T, const double x)
+{
+    // OGS 5 viscosity model 26
+
+    const double M0 = ProcessLib::TES::M_N2;
+    const double M1 = ProcessLib::TES::M_H2O;
+
+    // reactive component
+    const double x0 =
+        M0 * x / (M0 * x + M1 * (1.0 - x));  // mass in mole fraction
+    const double V0 = FluidViscosityH2O::get(M1 * p / (GAS_CONST * T), T);
+    // inert component
+    const double x1 = 1.0 - x0;
+    const double V1 = FluidViscosityN2::get(M0 * p / (GAS_CONST * T), T);
+
+    const double M0_over_M1(M1 / M0);  // reactive over inert
+    const double V0_over_V1(V0 / V1);
+
+    const double phi_12 =
+        mypow<2>(1.0 +
+                 std::sqrt(V0_over_V1) * std::pow(1.0 / M0_over_M1, 0.25)) /
+        std::sqrt(8.0 * (1.0 + M0_over_M1));
+
+    const double phi_21 = phi_12 * M0_over_M1 / V0_over_V1;
+
+    return V0 * x0 / (x0 + x1 * phi_12) + V1 * x1 / (x1 + x0 * phi_21);
+}
+
+struct FluidHeatConductivityN2
+{
+    static double get(double rho, double T)
+    {
+        const double X1 = 0.95185202;
+        const double X2 = 1.0205422;
+
+        const double rho_c = 314;  // [kg/m3]
+        const double M = 28.013;
+        const double k = 1.38062e-23;
+        const double eps = 138.08483e-23;
+        const double N_A = 6.02213E26;
+        const double R = 8.31434;
+        // const double R = GAS_CONST;
+        const double CCF = 4.173;  // mW/m/K
+
+        const double c1 = 0.3125;
+        const double c2 = 2.0442e-49;
+        const double sigma = 0.36502496e-09;
+
+        rho /= rho_c;
+
+        // dilute heat conductivity
+        const double sum1 = loop1_term<0>(T) + loop1_term<1>(T) +
+                            loop1_term<2>(T) + loop1_term<3>(T) +
+                            loop1_term<4>(T) + loop1_term<5>(T) +
+                            loop1_term<6>(T);
+        const double temp(std::exp((f[8] / T)) - 1);
+        const double c_v0 =
+            R *
+            (sum1 + ((f[7] * (f[8] / T) * (f[8] / T) * (std::exp((f[8] / T)))) /
+                         (temp * temp) -
+                     1));
+
+        double cvint;
+        cvint = c_v0 * 1000 / N_A;
+
+        // dilute gas viscosity
+        const double log_T_star = std::log(T * k / eps);
+
+        const double Omega =
+            std::exp(loop2_term<0>(log_T_star) + loop2_term<1>(log_T_star) +
+                     loop2_term<2>(log_T_star) + loop2_term<3>(log_T_star) +
+                     loop2_term<4>(log_T_star));
+
+        // eta in [Pa*s]
+        const double eta_0 =
+            1e6 * (c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega));
+
+        const double F = eta_0 * k * N_A / (M * 1000);
+
+        const double lambda_tr = 2.5 * (1.5 - X1);
+        const double lambda_in = X2 * (cvint / k + X1);
+
+        const double lambda_0 = F * (lambda_tr + lambda_in);
+
+        const double sum2 = loop3_term<0>(rho) + loop3_term<1>(rho) +
+                            loop3_term<2>(rho) + loop3_term<3>(rho);
+        const double lambda_r = sum2 * CCF;
+
+        return (lambda_0 + lambda_r) / 1000;  // lambda in [W/m/K]
+    }
+
+private:
+    template <int i>
+    static double loop1_term(const double T)
+    {
+        return f[i] * mypow<i - 3>(T);
+    }
+
+    template <int i>
+    static double loop2_term(const double log_T_star)
+    {
+        return A[i] * mypow<i>(log_T_star);
+    }
+
+    template <int i>
+    static double loop3_term(const double rho)
+    {
+        return C[i] * mypow<i + 1>(rho);
+    }
+
+    const static double A[5];
+    const static double f[9];
+    const static double C[4];
+};
+
+struct FluidHeatConductivityH2O
+{
+    static double get(double rho, double T)
+    {
+        double S, Q;
+        double b[3], B[2], d[4], C[6];
+
+        T /= 647.096;
+        rho /= 317.11;
+
+        b[0] = -0.397070;
+        b[1] = 0.400302;
+        b[2] = 1.060000;
+
+        B[0] = -0.171587;
+        B[1] = 2.392190;
+
+        d[0] = 0.0701309;
+        d[1] = 0.0118520;
+        d[2] = 0.00169937;
+        d[3] = -1.0200;
+
+        C[0] = 0.642857;
+        C[1] = -4.11717;
+        C[2] = -6.17937;
+        C[3] = 0.00308976;
+        C[4] = 0.0822994;
+        C[5] = 10.0932;
+
+        const double sum1 = loop_term<0>(T) + loop_term<1>(T) +
+                            loop_term<2>(T) + loop_term<3>(T);
+
+        const double lambda_0 = std::sqrt(T) * sum1;
+        const double lambda_1 =
+            b[0] + b[1] * rho +
+            b[2] * std::exp(B[0] * (rho + B[1]) * (rho + B[1]));
+
+        const double dT = fabs(T - 1) + C[3];
+        const double dT_pow_3_5 = std::pow(dT, 3. / 5.);
+        Q = 2 + (C[4] / dT_pow_3_5);
+
+        if (T >= 1)
+            S = 1 / dT;
+        else
+            S = C[5] / dT_pow_3_5;
+
+        const double rho_pow_9_5 = std::pow(rho, 9. / 5.);
+        const double rho_pow_Q = std::pow(rho, Q);
+        const double T_pow_3_2 = T * std::sqrt(T);
+        const double lambda_2 =
+            (d[0] / mypow<10>(T) + d[1]) * rho_pow_9_5 *
+                std::exp(C[0] * (1 - rho * rho_pow_9_5)) +
+            d[2] * S * rho_pow_Q *
+                std::exp((Q / (1. + Q)) * (1 - rho * rho_pow_Q)) +
+            d[3] * std::exp(C[1] * T_pow_3_2 + C[2] / mypow<5>(rho));
+
+        return lambda_0 + lambda_1 + lambda_2;  // lambda in [W/m/K]
+    }
+
+private:
+    template <unsigned i>
+    static double loop_term(const double T)
+    {
+        return a[i] * mypow<i>(T);
+    }
+
+    static const double a[4];
+};
+
+inline double fluid_heat_conductivity(const double p,
+                                      const double T,
+                                      const double x)
+{
+    // OGS 5 fluid heat conductivity model 11
+
+    const double M0 = ProcessLib::TES::M_N2;
+    const double M1 = ProcessLib::TES::M_H2O;
+
+    // TODO [CL] max() is redundant if the fraction is guaranteed to be between
+    // 0 and 1.
+    // reactive component
+    const double x0 = std::max(M0 * x / (M0 * x + M1 * (1.0 - x)),
+                               0.);  // convert mass to mole fraction
+    const double k0 =
+        FluidHeatConductivityH2O::get(M1 * p / (GAS_CONST * T), T);
+    // inert component
+    const double x1 = 1.0 - x0;
+    const double k1 = FluidHeatConductivityN2::get(M0 * p / (GAS_CONST * T), T);
+
+    const double M1_over_M2 = M1 / M0;  // reactive over inert
+    const double V1_over_V2 =
+        FluidViscosityH2O::get(M1 * p / (GAS_CONST * T), T) /
+        FluidViscosityN2::get(M0 * p / (GAS_CONST * T), T);
+    const double L1_over_L2 = V1_over_V2 / M1_over_M2;
+
+    const double M12_pow_mquarter = std::pow(M1_over_M2, -0.25);
+    const double phi_12 = (1.0 + std::sqrt(L1_over_L2) * M12_pow_mquarter) *
+                          (1.0 + std::sqrt(V1_over_V2) * M12_pow_mquarter) /
+                          std::sqrt(8.0 * (1.0 + M1_over_M2));
+    const double phi_21 = phi_12 * M1_over_M2 / V1_over_V2;
+
+    return k0 * x0 / (x0 + x1 * phi_12) + k1 * x1 / (x1 + x0 * phi_21);
+}
+
+}  // TES
+}  // ProcessLib
+
+#endif  // PROCESSLIB_TES_OGS5MATERIALMODELS

ProcessLib/TES/TESProcess.cpp

+/**
+ * \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 "TESProcess.h"
+
+// TODO Copied from VectorMatrixAssembler. Could be provided by the DOF table.
+inline AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices
+getRowColumnIndices_(std::size_t const id,
+                     AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                     std::vector<GlobalIndexType>& indices)
+{
+    assert(dof_table.size() > id);
+    indices.clear();
+
+    // Local matrices and vectors will always be ordered by component,
+    // no matter what the order of the global matrix is.
+    for (unsigned c = 0; c < dof_table.getNumComponents(); ++c)
+    {
+        auto const& idcs = dof_table(id, c).rows;
+        indices.reserve(indices.size() + idcs.size());
+        indices.insert(indices.end(), idcs.begin(), idcs.end());
+    }
+
+    return AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
+                                                                 indices);
+}
+
+template <typename GlobalVector>
+void getVectorValues(
+    GlobalVector const& x,
+    AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& r_c_indices,
+    std::vector<double>& local_x)
+{
+    auto const& indices = r_c_indices.rows;
+    local_x.clear();
+    local_x.reserve(indices.size());
+
+    for (auto i : indices)
+    {
+        local_x.emplace_back(x.get(i));
+    }
+}
+
+// TODO that essentially duplicates code which is also present in ProcessOutput.
+template <typename GlobalVector>
+double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
+                     AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                     std::size_t const node_id,
+                     std::size_t const global_component_id)
+{
+    MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
+                              node_id};
+
+    auto const index = dof_table.getLocalIndex(
+        l, global_component_id, x.getRangeBegin(), x.getRangeEnd());
+
+    return x.get(index);
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalSetup>
+TESProcess<GlobalSetup>::TESProcess(
+    MeshLib::Mesh& mesh,
+    typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+    std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+        time_discretization,
+    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+    ProcessOutput<GlobalVector>&& process_output,
+    const BaseLib::ConfigTree& config)
+    : Process<GlobalSetup>(
+          mesh, nonlinear_solver, std::move(time_discretization),
+          std::move(process_variables), std::move(secondary_variables),
+          std::move(process_output))
+{
+    DBUG("Create TESProcess.");
+
+    // physical parameters for local assembly
+    {
+        std::vector<std::pair<std::string, double*>> params{
+            {"fluid_specific_heat_source",
+             &_assembly_params.fluid_specific_heat_source},
+            {"fluid_specific_isobaric_heat_capacity", &_assembly_params.cpG},
+            {"solid_specific_heat_source",
+             &_assembly_params.solid_specific_heat_source},
+            {"solid_heat_conductivity", &_assembly_params.solid_heat_cond},
+            {"solid_specific_isobaric_heat_capacity", &_assembly_params.cpS},
+            {"tortuosity", &_assembly_params.tortuosity},
+            {"diffusion_coefficient",
+             &_assembly_params.diffusion_coefficient_component},
+            {"porosity", &_assembly_params.poro},
+            {"solid_density_dry", &_assembly_params.rho_SR_dry},
+            {"solid_density_initial", &_assembly_params.initial_solid_density}};
+
+        for (auto const& p : params)
+        {
+            if (auto const par = config.getConfParamOptional<double>(p.first))
+            {
+                DBUG("setting parameter `%s' to value `%g'", p.first.c_str(),
+                     *par);
+                *p.second = *par;
+            }
+        }
+    }
+
+    // characteristic values of primary variables
+    {
+        std::vector<std::pair<std::string, Trafo*>> const params{
+            {"characteristic_pressure", &_assembly_params.trafo_p},
+            {"characteristic_temperature", &_assembly_params.trafo_T},
+            {"characteristic_vapour_mass_fraction", &_assembly_params.trafo_x}};
+
+        for (auto const& p : params)
+        {
+            if (auto const par = config.getConfParamOptional<double>(p.first))
+            {
+                INFO("setting parameter `%s' to value `%g'", p.first.c_str(),
+                     *par);
+                *p.second = Trafo{*par};
+            }
+        }
+    }
+
+    // permeability
+    if (auto par =
+            config.getConfParamOptional<double>("solid_hydraulic_permeability"))
+    {
+        DBUG(
+            "setting parameter `solid_hydraulic_permeability' to isotropic "
+            "value `%g'",
+            *par);
+        const auto dim = mesh.getDimension();
+        _assembly_params.solid_perm_tensor =
+            Eigen::MatrixXd::Identity(dim, dim) * (*par);
+    }
+
+    // reactive system
+    _assembly_params.react_sys = Adsorption::AdsorptionReaction::newInstance(
+        config.getConfSubtree("reactive_system"));
+
+    // debug output
+    if (auto const param =
+            config.getConfParamOptional<bool>("output_element_matrices"))
+    {
+        DBUG("output_element_matrices: %s", (*param) ? "true" : "false");
+
+        _assembly_params.output_element_matrices = *param;
+    }
+
+    // TODO somewhere else
+    /*
+    if (auto const param =
+    config.getConfParamOptional<bool>("output_global_matrix"))
+    {
+        DBUG("output_global_matrix: %s", (*param) ? "true" : "false");
+
+        this->_process_output.output_global_matrix = *param;
+    }
+    */
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::initializeConcreteProcess(
+    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh, unsigned const integration_order)
+{
+    DBUG("Create global assembler.");
+    _global_assembler.reset(new GlobalAssembler(dof_table));
+
+    ProcessLib::createLocalAssemblers<GlobalSetup, TESLocalAssembler>(
+        mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
+        _local_assemblers, _assembly_params);
+
+    // TODO move the two data members somewhere else.
+    // for extrapolation of secondary variables
+    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
+        all_mesh_subsets_single_component;
+    all_mesh_subsets_single_component.emplace_back(
+        new MeshLib::MeshSubsets(this->_mesh_subset_all_nodes.get()));
+    _local_to_global_index_map_single_component.reset(
+        new AssemblerLib::LocalToGlobalIndexMap(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            AssemblerLib::ComponentOrder::BY_LOCATION));
+
+    _extrapolator.reset(
+        new ExtrapolatorImplementation(MathLib::MatrixSpecifications(
+            0u, 0u, nullptr, _local_to_global_index_map_single_component.get(),
+            &mesh)));
+
+    // secondary variables
+    auto add2nd = [&](std::string const& var_name, unsigned const n_components,
+                      SecondaryVariableFunctions<GlobalVector>&& fcts) {
+        this->_secondary_variables.addSecondaryVariable(var_name, n_components,
+                                                        std::move(fcts));
+    };
+    auto makeEx =
+        [&](TESIntPtVariables var) -> SecondaryVariableFunctions<GlobalVector> {
+        return ProcessLib::makeExtrapolator(var, *_extrapolator,
+                                            _local_assemblers);
+    };
+
+    add2nd("solid_density", 1, makeEx(TESIntPtVariables::SOLID_DENSITY));
+    add2nd("reaction_rate", 1, makeEx(TESIntPtVariables::REACTION_RATE));
+    add2nd("velocity_x", 1, makeEx(TESIntPtVariables::VELOCITY_X));
+    if (mesh.getDimension() >= 2)
+        add2nd("velocity_y", 1, makeEx(TESIntPtVariables::VELOCITY_Y));
+    if (mesh.getDimension() >= 3)
+        add2nd("velocity_z", 1, makeEx(TESIntPtVariables::VELOCITY_Z));
+
+    add2nd("loading", 1, makeEx(TESIntPtVariables::LOADING));
+    add2nd("reaction_damping_factor", 1,
+           makeEx(TESIntPtVariables::REACTION_DAMPING_FACTOR));
+
+    namespace PH = std::placeholders;
+    using Self = TESProcess<GlobalSetup>;
+
+    add2nd("vapour_partial_pressure", 1,
+           {std::bind(&Self::computeVapourPartialPressure, this, PH::_1, PH::_2,
+                      PH::_3),
+            nullptr});
+    add2nd("relative_humidity", 1, {std::bind(&Self::computeRelativeHumidity,
+                                              this, PH::_1, PH::_2, PH::_3),
+                                    nullptr});
+    add2nd("equilibrium_loading", 1,
+           {std::bind(&Self::computeEquilibriumLoading, this, PH::_1, PH::_2,
+                      PH::_3),
+            nullptr});
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::assembleConcreteProcess(const double t,
+                                                      GlobalVector const& x,
+                                                      GlobalMatrix& M,
+                                                      GlobalMatrix& K,
+                                                      GlobalVector& b)
+{
+    DBUG("Assemble TESProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalSetup::executeMemberDereferenced(*_global_assembler,
+                                           &GlobalAssembler::assemble,
+                                           _local_assemblers, t, x, M, K, b);
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::preTimestep(GlobalVector const& x, const double t,
+                                          const double delta_t)
+{
+    DBUG("new timestep");
+
+    _assembly_params.delta_t = delta_t;
+    _assembly_params.current_time = t;
+    ++_assembly_params.timestep;  // TODO remove that
+
+    _x_previous_timestep =
+        MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::preIteration(const unsigned iter,
+                                           GlobalVector const& /*x*/)
+{
+    _assembly_params.iteration_in_current_timestep = iter;
+    ++_assembly_params.total_iteration;
+    ++_assembly_params.number_of_try_of_iteration;
+}
+
+template <typename GlobalSetup>
+NumLib::IterationResult TESProcess<GlobalSetup>::postIteration(
+    GlobalVector const& x)
+{
+    if (this->_process_output.output_iteration_results)
+    {
+        DBUG("output results of iteration %li",
+             _assembly_params.total_iteration);
+        std::string fn =
+            "tes_iter_" + std::to_string(_assembly_params.total_iteration) +
+            +"_ts_" + std::to_string(_assembly_params.timestep) + "_" +
+            std::to_string(_assembly_params.iteration_in_current_timestep) +
+            "_" + std::to_string(_assembly_params.number_of_try_of_iteration) +
+            ".vtu";
+
+        this->output(fn, 0, x);
+    }
+
+    bool check_passed = true;
+
+    if (!Trafo::constrained)
+    {
+        // bounds checking only has to happen if the vapour mass fraction is
+        // non-logarithmic.
+
+        std::vector<GlobalIndexType> indices_cache;
+        std::vector<double> local_x_cache;
+        std::vector<double> local_x_prev_ts_cache;
+
+        auto check_variable_bounds = [&](std::size_t id,
+                                         LocalAssembler& loc_asm) {
+            auto const r_c_indices = getRowColumnIndices_(
+                id, *this->_local_to_global_index_map, indices_cache);
+            getVectorValues(x, r_c_indices, local_x_cache);
+            getVectorValues(*_x_previous_timestep, r_c_indices,
+                            local_x_prev_ts_cache);
+
+            if (!loc_asm.checkBounds(local_x_cache, local_x_prev_ts_cache))
+                check_passed = false;
+        };
+
+        GlobalSetup::executeDereferenced(check_variable_bounds,
+                                         _local_assemblers);
+    }
+
+    if (!check_passed)
+        return NumLib::IterationResult::REPEAT_ITERATION;
+
+    // TODO remove
+    DBUG("ts %lu iteration %lu (in current ts: %lu) try %u accepted",
+         _assembly_params.timestep, _assembly_params.total_iteration,
+         _assembly_params.iteration_in_current_timestep,
+         _assembly_params.number_of_try_of_iteration);
+
+    _assembly_params.number_of_try_of_iteration = 0;
+
+    return NumLib::IterationResult::SUCCESS;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeVapourPartialPressure(
+    typename TESProcess::GlobalVector const& x,
+    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<GlobalVector>& result_cache)
+{
+    assert(&dof_table == this->_local_to_global_index_map.get());
+
+    result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+         nullptr});
+
+    GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+    for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+    {
+        auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                     COMPONENT_ID_PRESSURE);
+        auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                        COMPONENT_ID_MASS_FRACTION);
+
+        auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+            x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+        // TODO Problems with PETSc? (local vs. global index)
+        result_cache->set(node_id, p * x_nV);
+    }
+
+    return *result_cache;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeRelativeHumidity(
+    typename TESProcess::GlobalVector const& x,
+    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<GlobalVector>& result_cache)
+{
+    assert(&dof_table == this->_local_to_global_index_map.get());
+
+    result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+         nullptr});
+
+    GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+    for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+    {
+        auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                     COMPONENT_ID_PRESSURE);
+        auto const T = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                     COMPONENT_ID_TEMPERATURE);
+        auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                        COMPONENT_ID_MASS_FRACTION);
+
+        auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+            x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+        auto const p_S =
+            Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(T);
+
+        // TODO Problems with PETSc? (local vs. global index)
+        result_cache->set(node_id, p * x_nV / p_S);
+    }
+
+    return *result_cache;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeEquilibriumLoading(
+    typename TESProcess::GlobalVector const& x,
+    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<GlobalVector>& result_cache)
+{
+    assert(&dof_table == this->_local_to_global_index_map.get());
+
+    result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+         nullptr});
+
+    GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+    for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+    {
+        auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                     COMPONENT_ID_PRESSURE);
+        auto const T = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                     COMPONENT_ID_TEMPERATURE);
+        auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+                                        COMPONENT_ID_MASS_FRACTION);
+
+        auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+            x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+        auto const p_V = p * x_nV;
+        auto const C_eq =
+            (p_V <= 0.0) ? 0.0
+                         : _assembly_params.react_sys->getEquilibriumLoading(
+                               p_V, T, _assembly_params.M_react);
+
+        // TODO Problems with PETSc? (local vs. global index)
+        result_cache->set(node_id, C_eq);
+    }
+
+    return *result_cache;
+}
+
+// Explicitly instantiate TESProcess for GlobalSetupType.
+template class TESProcess<GlobalSetupType>;
+
+}  // namespace TES
+
+}  // namespace ProcessLib

ProcessLib/TES/TESProcess.h

+/**
+ * \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 PROCESS_LIB_TESPROCESS_H_
+#define PROCESS_LIB_TESPROCESS_H_
+
+#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+
+#include "TESAssemblyParams.h"
+#include "TESLocalAssembler.h"
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalSetup>
+class TESProcess final : public Process<GlobalSetup>
+{
+    using BP = Process<GlobalSetup>;  //!< "Base Process"
+
+public:
+    using GlobalVector = typename GlobalSetup::VectorType;
+    using GlobalMatrix = typename GlobalSetup::MatrixType;
+
+    TESProcess(
+        MeshLib::Mesh& mesh,
+        typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+        std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+            time_discretization,
+        std::vector<std::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+        ProcessOutput<GlobalVector>&& process_output,
+        BaseLib::ConfigTree const& config);
+
+    void preTimestep(GlobalVector const& x, const double t,
+                     const double delta_t) override;
+    void preIteration(const unsigned iter, GlobalVector const& x) override;
+    NumLib::IterationResult postIteration(GlobalVector const& x) override;
+
+    bool isLinear() const override { return false; }
+private:
+    using LocalAssembler =
+        TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
+
+    using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
+        GlobalMatrix, GlobalVector, LocalAssembler,
+        NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
+
+    using ExtrapolatorInterface =
+        NumLib::Extrapolator<GlobalVector, TESIntPtVariables, LocalAssembler>;
+    using ExtrapolatorImplementation =
+        NumLib::LocalLinearLeastSquaresExtrapolator<
+            GlobalVector, TESIntPtVariables, LocalAssembler>;
+
+    void initializeConcreteProcess(
+        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh, unsigned const integration_order) override;
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
+
+    GlobalVector const& computeVapourPartialPressure(
+        GlobalVector const& x,
+        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        std::unique_ptr<GlobalVector>& result_cache);
+
+    GlobalVector const& computeRelativeHumidity(
+        GlobalVector const& x,
+        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        std::unique_ptr<GlobalVector>& result_cache);
+
+    GlobalVector const& computeEquilibriumLoading(
+        GlobalVector const& x,
+        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        std::unique_ptr<GlobalVector>& result_cache);
+
+    std::unique_ptr<GlobalAssembler> _global_assembler;
+    std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
+
+    AssemblyParams _assembly_params;
+
+    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_single_component;
+
+    std::unique_ptr<ExtrapolatorInterface> _extrapolator;
+
+    // used for checkBounds()
+    std::unique_ptr<GlobalVector> _x_previous_timestep;
+};
+
+template <typename GlobalSetup>
+std::unique_ptr<TESProcess<GlobalSetup>> createTESProcess(
+    MeshLib::Mesh& mesh,
+    typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+    std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+        time_discretization,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& /*parameters*/,
+    BaseLib::ConfigTree const& config)
+{
+    config.checkConfParam("type", "TES");
+
+    DBUG("Create TESProcess.");
+
+    auto process_variables = findProcessVariables(
+        variables, config,
+        {"fluid_pressure", "temperature", "vapour_mass_fraction"});
+
+    SecondaryVariableCollection<typename GlobalSetup::VectorType>
+        secondary_variables{
+            config.getConfSubtreeOptional("secondary_variables"),
+            {"solid_density", "reaction_rate", "velocity_x", "velocity_y",
+             "velocity_z", "loading", "reaction_damping_factor",
+             "vapour_partial_pressure", "relative_humidity",
+             "equilibrium_loading"}};
+
+    ProcessOutput<typename GlobalSetup::VectorType> process_output{
+        config.getConfSubtree("output"), process_variables,
+        secondary_variables};
+
+    return std::unique_ptr<TESProcess<GlobalSetup>>{new TESProcess<GlobalSetup>{
+        mesh, nonlinear_solver, std::move(time_discretization),
+        std::move(process_variables), std::move(secondary_variables),
+        std::move(process_output), config}};
+}
+
+}  // namespace TES
+
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_TESPROCESS_H_

ProcessLib/TES/TESReactionAdaptor.cpp

+/**
+ * \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 "MathLib/Nonlinear/Root1D.h"
+#include "MathLib/ODE/ODESolverBuilder.h"
+
+#include "MaterialsLib/Adsorption/Adsorption.h"
+#include "MaterialsLib/Adsorption/ReactionInert.h"
+#include "MaterialsLib/Adsorption/ReactionSinusoidal.h"
+
+#include "TESLocalAssemblerInner.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+std::unique_ptr<TESFEMReactionAdaptor> TESFEMReactionAdaptor::newInstance(
+    TESLocalAssemblerData const& data)
+{
+    auto const* ads = data.ap.react_sys.get();
+    if (dynamic_cast<Adsorption::AdsorptionReaction const*>(ads) != nullptr)
+    {
+        return std::unique_ptr<TESFEMReactionAdaptor>(
+            new TESFEMReactionAdaptorAdsorption(data));
+    }
+    else if (dynamic_cast<Adsorption::ReactionInert const*>(ads) != nullptr)
+    {
+        return std::unique_ptr<TESFEMReactionAdaptor>(
+            new TESFEMReactionAdaptorInert(data));
+    }
+    else if (dynamic_cast<Adsorption::ReactionSinusoidal const*>(ads) !=
+             nullptr)
+    {
+        return std::unique_ptr<TESFEMReactionAdaptor>(
+            new TESFEMReactionAdaptorSinusoidal(data));
+    }
+    else if (dynamic_cast<Adsorption::ReactionCaOH2 const*>(ads) != nullptr)
+    {
+        return std::unique_ptr<TESFEMReactionAdaptor>(
+            new TESFEMReactionAdaptorCaOH2(data));
+    }
+
+    ERR("No suitable TESFEMReactionAdaptor found. Aborting.");
+    std::abort();
+    return std::unique_ptr<TESFEMReactionAdaptor>(nullptr);
+}
+
+TESFEMReactionAdaptorAdsorption::TESFEMReactionAdaptorAdsorption(
+    TESLocalAssemblerData const& data)
+    // caution fragile: this relies in this constructor b eing called __after__
+    // data.solid_density has been properly set up!
+    : _bounds_violation(data.solid_density.size(), false),
+      _d(data)
+{
+    assert(dynamic_cast<Adsorption::AdsorptionReaction const*>(
+               data.ap.react_sys.get()) != nullptr &&
+           "Reactive system has wrong type.");
+    assert(_bounds_violation.size() != 0);
+}
+
+ReactionRate
+TESFEMReactionAdaptorAdsorption::initReaction_slowDownUndershootStrategy(
+    const unsigned int_pt)
+{
+    assert(_d.ap.number_of_try_of_iteration <= 20);
+
+    const double loading = Adsorption::AdsorptionReaction::getLoading(
+        _d.solid_density_prev_ts[int_pt], _d.ap.rho_SR_dry);
+
+    double react_rate_R =
+        _d.ap.react_sys->getReactionRate(_d.p_V, _d.T, _d.ap.M_react, loading) *
+        _d.ap.rho_SR_dry;
+
+    // set reaction rate based on current damping factor
+    react_rate_R = (_reaction_damping_factor > 1e-3)
+                       ? _reaction_damping_factor * react_rate_R
+                       : 0.0;
+
+    if (_d.p_V <
+            0.01 * Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(
+                       _d.T) &&
+        react_rate_R > 0.0)
+    {
+        react_rate_R = 0.0;
+    }
+    else if (_d.p_V < 100.0 ||
+             _d.p_V < 0.05 * Adsorption::AdsorptionReaction::
+                                 getEquilibriumVapourPressure(_d.T))
+    {
+        // use equilibrium reaction for dry regime
+
+        // in the case of zeroth try in zeroth iteration: _p_V and loading are
+        // the values
+        // at the end of the previous timestep
+
+        const double pV_eq = estimateAdsorptionEquilibrium(_d.p_V, loading);
+        // TODO [CL]: it would be more correct to subtract pV from the previous
+        // timestep here
+        const double delta_pV = pV_eq - _d.p_V;
+        const double delta_rhoV = delta_pV * _d.ap.M_react /
+                                  Adsorption::GAS_CONST / _d.T * _d.ap.poro;
+        const double delta_rhoSR = delta_rhoV / (_d.ap.poro - 1.0);
+        double react_rate_R2 = delta_rhoSR / _d.ap.delta_t;
+
+        if (_bounds_violation[int_pt])
+        {
+            react_rate_R2 *= 0.5;
+        }
+
+        // 0th try: make sure reaction is not slower than allowed by local
+        // estimation
+        // nth try: make sure reaction is not made faster by local estimation
+        if ((_d.ap.number_of_try_of_iteration == 1 &&
+             std::abs(react_rate_R2) > std::abs(react_rate_R)) ||
+            (_d.ap.number_of_try_of_iteration > 1 &&
+             std::abs(react_rate_R2) < std::abs(react_rate_R)))
+        {
+            react_rate_R = react_rate_R2;
+        }
+    }
+
+    // smooth out readjustment of reaction rate
+    if (_d.ap.iteration_in_current_timestep > 4)
+    {
+        if (_d.ap.iteration_in_current_timestep <= 9)
+        {
+            // update reaction rate for for five iterations
+            const auto N = _d.ap.iteration_in_current_timestep - 4;
+
+            // take average s.t. does not oscillate so much
+            react_rate_R = 1.0 / (1.0 + N) *
+                           (N * _d.reaction_rate[int_pt] + 1.0 * react_rate_R);
+        }
+        else
+        {
+            // afterwards no update anymore
+            react_rate_R = _d.reaction_rate[int_pt];
+        }
+    }
+
+    if (_d.ap.number_of_try_of_iteration > 1)
+    {
+        // assert that within tries reaction does not get faster
+        // (e.g. due to switch equilibrium reaction <--> kinetic reaction)
+
+        // factor of 0.9*N: in fact, even slow down reaction over tries
+        const double r = std::pow(0.9, _d.ap.number_of_try_of_iteration - 1) *
+                         _d.reaction_rate[int_pt];
+        if (std::abs(react_rate_R) > std::abs(r))
+        {
+            react_rate_R = r;
+        }
+    }
+
+    return {react_rate_R,
+            _d.solid_density_prev_ts[int_pt] + react_rate_R * _d.ap.delta_t};
+}
+
+double TESFEMReactionAdaptorAdsorption::estimateAdsorptionEquilibrium(
+    const double p_V0, const double C0) const
+{
+    auto f = [this, p_V0, C0](double pV) -> double {
+        // pV0 := _p_V
+        const double C_eq =
+            _d.ap.react_sys->getEquilibriumLoading(pV, _d.T, _d.ap.M_react);
+        return (pV - p_V0) * _d.ap.M_react / Adsorption::GAS_CONST / _d.T *
+                   _d.ap.poro +
+               (1.0 - _d.ap.poro) * (C_eq - C0) * _d.ap.rho_SR_dry;
+    };
+
+    // range where to search for roots of f
+    const double C_eq0 =
+        _d.ap.react_sys->getEquilibriumLoading(p_V0, _d.T, _d.ap.M_react);
+    const double limit =
+        (C_eq0 > C0)
+            ? 1e-8
+            : Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(
+                  _d.T);
+
+    // search for roots
+    auto rf = MathLib::Nonlinear::makeRegulaFalsi<MathLib::Nonlinear::Pegasus>(
+        f, p_V0, limit);
+    rf.step(3);
+
+    // set vapour pressure
+    return rf.getResult();
+}
+
+bool TESFEMReactionAdaptorAdsorption::checkBounds(
+    std::vector<double> const& local_x,
+    std::vector<double> const& local_x_prev_ts)
+{
+    double alpha = 1.0;
+
+    const double min_xmV = 1e-6;
+    const std::size_t nnodes = local_x.size() / NODAL_DOF;
+    const std::size_t xmV_offset = COMPONENT_ID_MASS_FRACTION * nnodes;
+
+    for (std::size_t i = 0; i < nnodes; ++i)
+    {
+        auto const xnew = local_x[xmV_offset + i];
+        auto const xold = local_x_prev_ts[xmV_offset + i];
+
+        if (xnew < min_xmV)
+        {
+            const auto a = xold / (xold - xnew);
+            alpha = std::min(alpha, a);
+            _bounds_violation[i] = true;
+        }
+        else if (xnew > 1.0)
+        {
+            const auto a = xold / (xnew - xold);
+            alpha = std::min(alpha, a);
+            _bounds_violation[i] = true;
+        }
+        else
+        {
+            _bounds_violation[i] = false;
+        }
+    }
+
+    assert(alpha > 0.0);
+
+    if (alpha != 1.0)
+    {
+        if (alpha > 0.5)
+            alpha = 0.5;
+        if (alpha < 0.05)
+            alpha = 0.05;
+
+        if (_d.ap.number_of_try_of_iteration <= 3)
+        {
+            _reaction_damping_factor *= sqrt(alpha);
+        }
+        else
+        {
+            _reaction_damping_factor *= alpha;
+        }
+    }
+
+    return alpha == 1.0;
+}
+
+void TESFEMReactionAdaptorAdsorption::preZerothTryAssemble()
+{
+    if (_reaction_damping_factor < 1e-3)
+        _reaction_damping_factor = 1e-3;
+
+    _reaction_damping_factor = std::min(std::sqrt(_reaction_damping_factor),
+                                        10.0 * _reaction_damping_factor);
+}
+
+TESFEMReactionAdaptorInert::TESFEMReactionAdaptorInert(
+    TESLocalAssemblerData const& data)
+    : _d(data)
+{
+}
+
+ReactionRate TESFEMReactionAdaptorInert::initReaction(const unsigned int_pt)
+{
+    return {0.0, _d.solid_density_prev_ts[int_pt]};
+    // _d._qR = 0.0;
+    // _d._reaction_rate[int_pt] = 0.0;
+}
+
+TESFEMReactionAdaptorSinusoidal::TESFEMReactionAdaptorSinusoidal(
+    TESLocalAssemblerData const& data)
+    : _d(data)
+{
+    assert(dynamic_cast<Adsorption::ReactionSinusoidal const*>(
+               data.ap.react_sys.get()) != nullptr &&
+           "Reactive system has wrong type.");
+}
+
+ReactionRate TESFEMReactionAdaptorSinusoidal::initReaction(
+    const unsigned /*int_pt*/)
+{
+    const double t = _d.ap.current_time;
+
+    // Cf. OGS5
+    const double rhoSR0 = 1.0;
+    const double rhoTil = 0.1;
+    const double omega = 2.0 * 3.1416;
+    const double poro = _d.ap.poro;
+
+    return {rhoTil * omega * cos(omega * t) / (1.0 - poro),
+            rhoSR0 + rhoTil * std::sin(omega * t) / (1.0 - poro)};
+}
+
+TESFEMReactionAdaptorCaOH2::TESFEMReactionAdaptorCaOH2(
+    TESLocalAssemblerData const& data)
+    : _d(data),
+      _react(dynamic_cast<Adsorption::ReactionCaOH2&>(*data.ap.react_sys.get()))
+{
+    _ode_solver = MathLib::ODE::createODESolver<1>(_react.getOdeSolverConfig());
+    // TODO invalidate config
+
+    _ode_solver->setTolerance(1e-10, 1e-10);
+
+    auto f = [this](const double /*t*/,
+                    MathLib::ODE::MappedConstVector<1> const y,
+                    MathLib::ODE::MappedVector<1>
+                        ydot) -> bool {
+        ydot[0] = _react.getReactionRate(y[0]);
+        return true;
+    };
+
+    _ode_solver->setFunction(f, nullptr);
+}
+
+ReactionRate TESFEMReactionAdaptorCaOH2::initReaction(const unsigned int int_pt)
+{
+    // TODO if the first holds, the second also has to hold
+    if (_d.ap.iteration_in_current_timestep > 1 ||
+        _d.ap.number_of_try_of_iteration > 1)
+    {
+        return {_d.reaction_rate[int_pt], _d.solid_density[int_pt]};
+    }
+
+    // TODO: double check!
+    // const double xv_NR  = SolidProp->non_reactive_solid_volume_fraction;
+    // const double rho_NR = SolidProp->non_reactive_solid_density;
+    const double xv_NR = 0.0;
+    const double rho_NR = 0.0;
+
+    const double t0 = 0.0;
+    const double y0 =
+        (_d.solid_density_prev_ts[int_pt] - xv_NR * rho_NR) / (1.0 - xv_NR);
+
+    const double t_end = _d.ap.delta_t;
+
+    _react.updateParam(_d.T, _d.p, _d.vapour_mass_fraction,
+                       _d.solid_density_prev_ts[int_pt]);
+
+    _ode_solver->setIC(t0, {y0});
+    _ode_solver->preSolve();
+    _ode_solver->solve(t_end);
+
+    const double time_reached = _ode_solver->getTime();
+    (void)time_reached;
+    assert(std::abs(t_end - time_reached) <
+           std::numeric_limits<double>::epsilon());
+
+    auto const& y_new = _ode_solver->getSolution();
+    auto const& y_dot_new = _ode_solver->getYDot(t_end, y_new);
+
+    double rho_react;
+
+    // cut off when limits are reached
+    if (y_new[0] < _react.rho_low)
+        rho_react = _react.rho_low;
+    else if (y_new[0] > _react.rho_up)
+        rho_react = _react.rho_up;
+    else
+        rho_react = y_new[0];
+
+    return {y_dot_new[0] * (1.0 - xv_NR),
+            (1.0 - xv_NR) * rho_react + xv_NR * rho_NR};
+}
+
+}  // namespace TES
+}  // namespace ProcessLib